Iterative fourier reconstruction for laser surgery and other optical applications

ABSTRACT

Methods, systems and software for determine an optical surface model for an optical tissue system using Fourier transformation algorithms. A method of reconstructing optical tissues of an eye comprises transmitting an image through the optical tissues of the eye. The surface gradients from the transmitted image are measured across the optical tissues of the eye. A Fourier transform algorithm is applied to the surface gradients to reconstruct an optical surface model that corresponds to the optical tissues of the eye.

CROSS-REFERENCES TO RELATED APPLICATIONS

The present application is a continuation-in-part of U.S. patentapplication Ser. No. 10/601,048 filed Jun. 20, 2003 (Attorney Docket No.018158-021800US) and is related to U.S. Provisional Patent ApplicationNo. 60/550,514 filed Mar. 3, 2004 (Attorney Docket No. 018158-024300US),the entire disclosures of which are incorporated herein by reference forall purposes.

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY SPONSOREDRESEARCH OR DEVELOPMENT

NOT APPLICABLE

REFERENCE TO A “SEQUENCE LISTING,” A TABLE, OR A COMPUTER PROGRAMLISTING APPENDIX SUBMITTED ON A COMPACT DISK

NOT APPLICABLE

BACKGROUND OF THE INVENTION

The present invention generally relates to measuring optical errors ofoptical systems. More particularly, the invention relates to improvedmethods and systems for determining an optical surface model for anoptical tissue system of an eye, to improved methods and systems forreconstructing a wavefront surface/elevation map of optical tissues ofan eye, and to improved systems for calculating an ablation pattern.

Known laser eye surgery procedures generally employ an ultraviolet orinfrared laser to remove a microscopic layer of stromal tissue from thecornea of the eye. The laser typically removes a selected shape of thecorneal tissue, often to correct refractive errors of the eye.Ultraviolet laser ablation results in photodecomposition of the cornealtissue, but generally does not cause significant thermal damage toadjacent and underlying tissues of the eye. The irradiated molecules arebroken into smaller volatile fragments photochemically, directlybreaking the intermolecular bonds.

Laser ablation procedures can remove the targeted stroma of the corneato change the cornea's contour for varying purposes, such as forcorrecting myopia, hyperopia, astigmatism, and the like. Control overthe distribution of ablation energy across the cornea may be provided bya variety of systems and methods, including the use of ablatable masks,fixed and moveable apertures, controlled scanning systems, eye movementtracking mechanisms, and the like. In known systems, the laser beamoften comprises a series of discrete pulses of laser light energy, withthe total shape and amount of tissue removed being determined by theshape, size, location, and/or number of laser energy pulses impinging onthe cornea. A variety of algorithms may be used to calculate the patternof laser pulses used to reshape the cornea so as to correct a refractiveerror of the eye. Known systems make use of a variety of forms of lasersand/or laser energy to effect the correction, including infrared lasers,ultraviolet lasers, femtosecond lasers, wavelength multipliedsolid-state lasers, and the like. Alternative vision correctiontechniques make use of radial incisions in the cornea, intraocularlenses, removable corneal support structures, and the like.

Known corneal correction treatment methods have generally beensuccessful in correcting standard vision errors, such as myopia,hyperopia, astigmatism, and the like. However, as with all successes,still further improvements would be desirable. Toward that end,wavefront measurement systems are now available to accurately measurethe refractive characteristics of a particular patient's eye. Oneexemplary wavefront technology system is the VISX WaveScan® System,which uses a Hartmann-Shack wavefront lenslet array that can quantifyaberrations throughout the entire optical system of the patient's eye,including first- and second-order sphero-cylindrical errors, coma, andthird and fourth-order aberrations related to coma, astigmatism, andspherical aberrations.

Wavefront measurement of the eye may be used to create an ocularaberration map, a high order aberration map, or wavefront elevation mapthat permits assessment of aberrations throughout the optical pathway ofthe eye, e.g., both internal aberrations and aberrations on the cornealsurface. The aberration map may then be used to compute a customablation pattern for allowing a surgical laser system to correct thecomplex aberrations in and on the patient's eye. Known methods forcalculation of a customized ablation pattern using wavefront sensor datagenerally involve mathematically modeling an optical surface of the eyeusing expansion series techniques.

Reconstruction of the wavefront or optical path difference (OPD) of thehuman ocular aberrations can be beneficial for a variety of uses. Forexample, the wavefront map, the wavefront refraction, the point spreadfunction, and the treatment table can all depend on the reconstructedwavefront.

Known wavefront reconstruction can be categorized into two approaches:zonal reconstruction and modal reconstruction. Zonal reconstruction wasused in early adaptive optics systems. More recently, modalreconstruction has become popular because of the use of Zernikepolynomials. Coefficients of the Zernike polynomials can be derivedthrough known fitting techniques, and the refractive correctionprocedure can be determined using the shape of the optical surface ofthe eye, as indicated by the mathematical series expansion model.

The Zernike function method of surface reconstruction and its accuracyfor normal eyes have limits. For example, 6th order Zernike polynomialsmay not accurately represent an actual wavefront in all circumstances.The discrepancy may be most significant for eyes with a keratoconuscondition. Known Zernike polynomial modeling methods may also result inerrors or “noise” which can lead to a less than ideal refractivecorrection. Furthermore, the known surface modeling techniques aresomewhat indirect, and may lead to unnecessary errors in calculation, aswell as a lack of understanding of the physical correction to beperformed.

Therefore, in light of above, it would be desirable to provide improvedmethods and systems for mathematically modeling optical tissues of aneye.

BRIEF SUMMARY OF THE INVENTION

The present invention provides novel iterative Fourier transform methodsand systems that can account for missing, erroneous, or otherwiseinsufficient data points. The present invention also provides determinedgoals, or metrics that can be used to determine an optimum or reasonablenumber of iterations. What is more, the present invention providessystems, software, and methods for measuring errors and reconstructingwavefront elevation maps in an optical system using Fourier transformalgorithms.

In a first aspect, the present invention provides a method ofdetermining an optical surface model for an optical tissue system of aneye. The method can include a) inputting optical data from the opticaltissue system of the eye, the optical data comprising set of localgradients; b) establishing a first combined gradient field based on theset of local gradients; c) deriving a first reconstructed wavefront byapplying a Fourier transform to the first combined gradient field; d)providing a first revised gradient field based on the firstreconstructed wavefront; e) establishing a second combined gradientfield based on the first revised gradient field; f) deriving a secondreconstructed wavefront by applying the Fourier transform to the secondcombined gradient field; and g) determining the optical surface modelbased on the second reconstructed wavefront.

In some aspects, the optical data includes wavefront data. In someaspects, the wavefront data includes a set of local or surface gradientswithin an aperture. In some aspects, the aperture corresponds to a pupilof an eye.

In further aspects, the first combined gradient field includes a firstexterior gradient field and a measured gradient field, such that themeasured gradient field is disposed interior to the first exteriorgradient field, and the measured gradient field corresponds to the setof local gradients; and the second combined gradient field includes asecond exterior gradient field and a measured gradient field, such thatthe second exterior gradient field corresponds to the first revisedgradient field, and the measured gradient field is disposed interior tothe second exterior gradient field.

In related aspects, the method of determining an optical surface modelfor an optical tissue system of an eye also includes selecting at leasta portion of the second reconstructed wavefront to provide a finalreconstructed wavefront, and determining the optical surface model basedon the final reconstructed wavefront.

In related aspects, the method further comprises iteratively performingsteps (d) through (f) as noted above to derive an nth reconstructedwavefront by applying the Fourier transform to an nth combined gradientfield; selecting at least a portion of the nth reconstructed wavefrontto provide the final reconstructed wavefront; and determining theoptical surface model based on the final reconstructed wavefront.

In a related aspect, the present invention provides a method of mappinga wavefront error of an eye. The method includes determining an opticalsurface model as described above, and mapping the wavefront error of theeye based on the optical surface model. In a further related aspect, thepresent invention provides a method of computing a correction ablationpattern for an optical tissue system of an eye. The method includesdetermining an optical surface model as described above, and computingthe correction ablation pattern for the eye based on the optical surfacemodel.

In further aspects, the first combined gradient field includes a firstexterior gradient field and a measured gradient field, such that themeasured gradient field is disposed interior to the first exteriorgradient field, and the measured gradient field corresponds to the setof local gradients; and the second combined gradient field includes asecond exterior gradient field and a measured gradient field, such thatthe second exterior gradient field corresponds to the first revisedgradient field, and the measured gradient field is disposed interior tothe second exterior gradient field.

In related aspects, the method of determining an optical surface modelfor an optical tissue system of an eye also includes selecting at leasta portion of the second reconstructed wavefront to provide a finalreconstructed wavefront, and determining the optical surface model basedon the final reconstructed wavefront.

In a related aspect, the method can include iteratively performing steps(d) through (f) as described above to derive an nth reconstructedwavefront based on the application of the Fourier transform to an nthcombined gradient field; selecting at least a portion of the nthreconstructed wavefront to provide the final reconstructed wavefront;and determining the optical surface model based on the finalreconstructed wavefront.

In a related aspect, the present invention provides a method ofmodifying an optical tissue surface of an eye. The method can includecomputing a correction ablation pattern as described above; andmodifying the optical tissue surface according to the correctionablation pattern by applying a laser ablation to the eye.

In a second aspect, the present invention provides a system fordetermining an optical surface model for an optical tissue system of aneye. The system can include a) a light source for transmitting an imagethrough the optical tissue system; b) a sensor oriented for determininga set of local gradients for the optical tissue system by detecting thetransmitted image; c) a processor coupled to the sensor; and d) a memorycoupled to the processor, the memory configured to store a plurality ofcode modules for execution by the processor. The plurality of codemodules can include i) a module for establishing a first combinedgradient field based on the set of local gradients; ii) a module forderiving a first reconstructed wavefront by applying a Fourier transformto the first combined gradient field; iii) a module for providing afirst revised gradient field based on the first reconstructed wavefront;iv) a module for establishing a second combined gradient field based onthe first revised gradient field; v) a module for deriving a secondreconstructed wavefront by applying the Fourier transform to the secondcombined gradient field; and vi) a module for determining the opticalsurface model for the optical tissue system of the eye based on thesecond reconstructed wavefront. The modules described in this aspect, aswell as those described throughout the application, may include dataprocessing software and/or hardware, and may be integrated with otherdata processing structures.

In a related aspect, the invention provides a system for mapping awavefront error of an optical tissue system of an eye. The system caninclude a system as describe above for determining an optical surfacemodel for an optical tissue system of an eye, along with a module formapping a wavefront error of the eye based on the optical surface model.In a further related aspect, the present invention provides a system forcomputing a correction ablation pattern for an optical tissue system ofan eye. The system can include a system as described above fordetermining an optical surface model for an optical tissue system of aneye, along with a module for computing a correction ablation pattern forthe eye based on the optical surface model. In yet another aspect, theinvention provides a system for modifying an optical tissue surface ofan eye. The system can include a system as described above for computinga correction ablation pattern for an optical tissue system of an eye,along with a laser for modifying the optical tissue surface of the eyebased on the correction ablation pattern.

In a third aspect, the present invention provides a system fordetermining an optical surface model that corresponds to an opticaltissue system of an eye. The system can include a) means fortransmitting an image through the optical tissue system; b) means, in anoptical path from the image, for determining a set of local gradientsfor the optical tissue system based on the transmitted image; c) means,coupled to the local gradient determining means, for establishing afirst combined gradient field based on the set of local gradients; d)means, coupled to the first combined gradient field establishing means,for deriving a first reconstructed wavefront based on the application ofa Fourier transform to the first combined gradient field; e) means,coupled to the first reconstructed wavefront deriving means, forproviding a first revised gradient field based on the firstreconstructed wavefront; f) means, coupled to the first revised gradientfield providing means, for establishing a second combined gradient fieldbased on the first revised gradient field; g) means, coupled to thesecond combined gradient field establishing means, for deriving a secondreconstructed wavefront based on the application of the Fouriertransform to the second combined gradient field; and h) means, coupledto the second reconstructed wavefront deriving means, for determiningthe optical surface model for the optical tissue system of the eye basedon the second reconstructed wavefront.

In an fourth aspect, the present invention provides a computer programstored on a computer-readable storage medium. The computer program caninclude a) code for transmitting an image through an optical tissuesystem of an eye; b) code for determining a set of local gradients forthe optical tissue system of the eye based on the transmitted image; c)code for establishing a first combined gradient field based on the setof local gradients; d) code for deriving a first reconstructed wavefrontbased on the application of a Fourier transform to the first combinedgradient field; e) code for providing a first revised gradient fieldbased on the first reconstructed wavefront; f) code for establishing asecond combined gradient field based on the first revised gradientfield; g) code for deriving a second reconstructed wavefront based onthe application of the Fourier transform to the second combined gradientfield; and h) code for determining an optical surface model for theoptical tissue system of the eye based on the second reconstructedwavefront. In a related aspect, the computer program can also includecode for computing a correction ablation pattern based on the opticalsurface model. In a further related aspect, the computer program canalso include code for delivering a laser energy to the eye based on thecorrection ablation pattern.

In one aspect, the present invention provides a method of reconstructingoptical tissues of an eye. The method comprises transmitting an imagethrough the optical tissues of the eye. Surface gradients from thetransmitted image are measured across the optical tissues of the eye. AFourier transform algorithm is applied to the surface gradients toreconstruct a surface that corresponds to the optical tissues of theeye.

In some embodiments, the measurement of the surface gradient data iscarried out with a Hartmann-Shack sensor assembly. The image istransmitted by the optical tissues as a plurality of beamlets and thesurface gradients will be in the form of an array of gradients. Eachgradient corresponds to an associated portion of the optical tissues ofthe eye and each beamlet is transmitted through the optical tissuesaccording to the corresponding gradient. In such embodiments, themeasured surface will be in the form of a wavefront surface or wavefrontelevation map.

In one embodiment, the Fourier transformation algorithm may apply a fastFourier transform or a discrete Fourier decomposition and an inversediscrete Fourier transform. Some Fourier transform algorithms mayinclude a mean gradient field so as to remove a tilt from thereconstructed surface. Unlike Zernike polynomial reconstruction, whichis based on a unit circle, the Fourier transformation uses all of theavailable information in the reconstruction.

A computed correction ablation pattern based on the reconstructedoptical tissues of the eye, as indicated by the Fourier reconstructedsurface, may be calculated and aligned with the eye. The correctionablation pattern typically comprises a proposed change in elevations ofthe optical tissue so as to effect a desired change in opticalproperties of the eye. After the correction ablation pattern is derived,laser ablation may be used to modify the optical tissue surface.

In another aspect, the present invention provides a method for measuringoptical tissues of an eye. The method comprises transmitting an imagethrough the optical tissues. Local gradients across the optical tissuesare determined from the transmitted image. A wavefront error of the eyeis mapped by applying a Fourier transform algorithm to the surfacegradients across the optical tissues of the eye.

Measurement of the local gradients may be carried out by aHartmann-Shack sensor assembly. A mean gradient field may be added tothe wavefront error to correct for tilt. After the wavefront error ismapped, a laser ablation treatment table may be created based on themapped wavefront error of the optical tissues of the eye and the opticaltissue surface may be modified according to the correction ablationpattern by laser ablation.

In another aspect, the present invention further provides a system formeasuring a wavefront error of optical tissue. The system comprises amemory coupled to a processor. The memory may be configured to store aplurality of code modules for execution by the processor. The pluralityof code modules comprise a module for transmitting an image through theoptical tissues, a module for determining local gradients across theoptical tissues from the transmitted image, and a module for mapping awavefront error of the eye by applying a Fourier transform algorithm tothe surface gradients across the optical tissues of the eye.

The system may include an image source coupled to the processor fortransmitting a source image through the optical tissues of the eye. Theimage may be a point or small spot of light, or any other suitableimage. The system may also include a wavefront sensor system coupled tothe processor, such as a Hartmann-Shack sensor assembly.

The system may include one or more cameras to track the position of theoptical tissues. In such embodiments, the plurality of code modules alsoincludes a module for registering the wavefront error relative to animage of the optical tissues that was obtained by the camera(s).

In some embodiments, the system may include a module for calculating acustomized laser ablation program or ablation table based on thereconstructed surface of the optical tissues. A laser system may be incommunication with the above measurement system. The laser system mayinclude a laser that is programmable to deliver a laser energy to theoptical tissues according to the correction ablation pattern.

In a further aspect, the present invention provides a computer programstored on a computer-readable storage medium. The computer programcomprises code for transmitting an image through the optical tissues ofthe eye, code for measuring surface gradients from the transmitted imageacross the optical tissues of the eye, and code for mapping a wavefronterror of the eye by applying a Fourier transform algorithm to thesurface gradients across the optical tissues of the eye.

In some embodiments, the computer program may include code for computingan ablation pattern based on the optical tissues of the eye as indicatedby the Fourier reconstructed surface, code for controlling the deliveryof laser energy to the optical tissues according to the correctionablation pattern, and/or code for aligning the mapped wavefront errorwith an image of the optical tissues of the eye.

The methods and apparatuses of the present invention may be provided inone or more kits for such use. For example, the kits may comprise asystem for determining an optical surface model that corresponds to anoptical tissue system of an eye. Optionally, such kits may furtherinclude any of the other system components described in relation to thepresent invention and any other materials or items relevant to thepresent invention. The instructions for use can set forth any of themethods as described above. It is further understood that systemsaccording to the present invention may be configured to carry out any ofthe method steps described above.

These and other aspects will be apparent in the remainder of thefigures, description and claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a laser ablation system according to an embodiment ofthe present invention.

FIG. 2 illustrates a simplified computer system according to anembodiment of the present invention.

FIG. 3 illustrates a wavefront measurement system according to anembodiment of the present invention.

FIG. 3A illustrates another wavefront measurement system according toanother embodiment of the present invention.

FIG. 4 schematically illustrates a simplified set of modules that carryout one method of the present invention.

FIG. 5 is a flow chart that schematically illustrates a method of usinga Fourier transform algorithm to determine a corneal ablation treatmentprogram.

FIG. 6 schematically illustrates a comparison of a direct integrationreconstruction, a 6th order Zernike polynomial reconstruction, a 10thorder Zernike polynomial reconstruction, and a Fourier transformreconstruction for a surface having a +2 ablation on a 6 mm pupil.

FIG. 7 schematically illustrates a comparison of a direct integrationreconstruction, a 6th order Zernike polynomial reconstruction, a 10thorder Zernike polynomial reconstruction, and a Fourier transformreconstruction for a presbyopia surface.

FIG. 8 schematically illustrates a comparison of a direct integrationreconstruction, a 6th order Zernike polynomial reconstruction, a 10thorder Zernike polynomial reconstruction, and a Fourier transformreconstruction for another presbyopia surface.

FIG. 9 illustrates a difference in a gradient field calculated from areconstructed wavefront from a Fourier transform reconstructionalgorithm (F Gradient), Zernike polynomial reconstruction algorithm (ZGradient), a direct integration reconstruction algorithm (D Gradient)and a directly measured wavefront.

FIG. 10 illustrates intensity plots of reconstructed wavefronts forFourier, 10th Order Zernike and Direct Integration reconstructions.

FIG. 11 illustrates intensity plots of reconstructed wavefronts forFourier, 6th Order Zernike and Direct Integration reconstructions.

FIG. 12 illustrates an algorithm flow chart according to one embodimentof the present invention.

FIG. 13 illustrates surface plots of wavefront reconstruction accordingto one embodiment of the present invention.

FIG. 14 illustrates surface plots of wavefront reconstruction accordingto one embodiment of the present invention.

FIG. 15 illustrates a comparison of wavefront maps with or withoutmissing data.

DETAILED DESCRIPTION OF THE INVENTION

The present invention provides systems, software, and methods that use ahigh speed and accurate Fourier or iterative Fourier transformationalgorithm to mathematically determine an optical surface model for anoptical tissue system of an eye or to otherwise mathematicallyreconstruct optical tissues of an eye.

The present invention is generally useful for enhancing the accuracy andefficacy of laser eye surgical procedures, such as photorefractivekeratectomy (PRK), phototherapeutic keratectomy (PTK), laser in situkeratomileusis (LASIK), and the like. The present invention can provideenhanced optical accuracy of refractive procedures by improving themethodology for measuring the optical errors of the eye and hencecalculate a more accurate refractive ablation program. In one particularembodiment, the present invention is related to therapeuticwavefront-based ablations of pathological eyes.

The present invention can be readily adapted for use with existing lasersystems, wavefront measurement systems, and other optical measurementdevices. By providing a more direct (and hence, less prone to noise andother error) methodology for measuring and correcting errors of anoptical system, the present invention may facilitate sculpting of thecornea so that treated eyes regularly exceed the normal 20/20 thresholdof desired vision. While the systems, software, and methods of thepresent invention are described primarily in the context of a laser eyesurgery system, it should be understood the present invention may beadapted for use in alternative eye treatment procedures and systems suchas spectacle lenses, intraocular lenses, contact lenses, corneal ringimplants, collagenous corneal tissue thermal remodeling, and the like.

Referring now to FIG. 1, a laser eye surgery system 10 of the presentinvention includes a laser 12 that produces a laser beam 14. Laser 12 isoptically coupled to laser delivery optics 16, which directs laser beam14 to an eye of patient P. A delivery optics support structure (notshown here for clarity) extends from a frame 18 supporting laser 12. Amicroscope 20 is mounted on the delivery optics support structure, themicroscope often being used to image a cornea of the eye.

Laser 12 generally comprises an excimer laser, ideally comprising anargon-fluorine laser producing pulses of laser light having a wavelengthof approximately 193 run. Laser 12 will preferably be designed toprovide a feedback stabilized fluence at the patient's eye, deliveredvia laser delivery optics 16. The present invention may also be usefulwith alternative sources of ultraviolet or infrared radiation,particularly those adapted to controllably ablate the corneal tissuewithout causing significant damage to adjacent and/or underlying tissuesof the eye. In alternate embodiments, the laser beam source employs asolid state laser source having a wavelength between 193 and 215 nm asdescribed in U.S. Pat. Nos. 5,520,679, and 5,144,630 to Lin and U.S.Pat. No. 5,742,626 to Mead, the full disclosures of which areincorporated herein by reference. In another embodiment, the lasersource is an infrared laser as described in U.S. Pat. Nos. 5,782,822 and6,090,102 to Telfair, the full disclosures of which are incorporatedherein by reference. Hence, although an excimer laser is theillustrative source of an ablating beam, other lasers may be used in thepresent invention.

Laser 12 and laser delivery optics 16 will generally direct laser beam14 to the eye of patient P under the direction of a computer system 22.Computer system 22 will often selectively adjust laser beam 14 to exposeportions of the cornea to the pulses of laser energy so as to effect apredetermined sculpting of the cornea and alter the refractivecharacteristics of the eye. In many embodiments, both laser 12 and thelaser delivery optical system 16 will be under control of computersystem 22 to effect the desired laser sculpting process, with thecomputer system effecting (and optionally modifying) the pattern oflaser pulses. The pattern of pulses may be summarized in machinereadable data of tangible media 29 in the form of a treatment table, andthe treatment table may be adjusted according to feedback input intocomputer system 22 from an automated image analysis system (or manuallyinput into the processor by a system operator) in response to real-timefeedback data provided from an ablation monitoring system feedbacksystem. The laser treatment system 10, and computer system 22 maycontinue and/or terminate a sculpting treatment in response to thefeedback, and may optionally also modify the planned sculpting based atleast in part on the feedback.

Additional components and subsystems may be included with laser system10, as should be understood by those of skill in the art. For example,spatial and/or temporal integrators may be included to control thedistribution of energy within the laser beam, as described in U.S. Pat.No. 5,646,791, the full disclosure of which is incorporated herein byreference. Ablation effluent evacuators/filters, aspirators, and otherancillary components of the laser surgery system are known in the art.Further details of suitable systems for performing a laser ablationprocedure can be found in commonly assigned U.S. Pat. Nos. 4,665,913,4,669,466, 4,732,148, 4,770,172, 4,773,414, 5,207,668, 5,108,388,5,219,343, 5,646,791 and 5,163,934, the complete disclosures of whichare incorporated herein by reference. Suitable systems also includecommercially available refractive laser systems such as thosemanufactured and/or sold by Alcon, Bausch & Lomb, Nidek, WaveLight,LaserSight, Schwind, Zeiss-Meditec, and the like.

FIG. 2 is a simplified block diagram of an exemplary computer system 22that may be used by the laser surgical system 10 of the presentinvention. Computer system 22 typically includes at least one processor52 which may communicate with a number of peripheral devices via a bussubsystem 54. These peripheral devices may include a storage subsystem56, comprising a memory subsystem 58 and a file storage subsystem 60,user interface input devices 62, user interface output devices 64, and anetwork interface subsystem 66. Network interface subsystem 66 providesan interface to outside networks 68 and/or other devices, such as thewavefront measurement system 30.

User interface input devices 62 may include a keyboard, pointing devicessuch as a mouse, trackball, touch pad, or graphics tablet, a scanner,foot pedals, a joystick, a touchscreen incorporated into the display,audio input devices such as voice recognition systems, microphones, andother types of input devices. User input devices 62 will often be usedto download a computer executable code from a tangible storage media 29embodying any of the methods of the present invention. In general, useof the term “input device” is intended to include a variety ofconventional and proprietary devices and ways to input information intocomputer system 22.

User interface output devices 64 may include a display subsystem, aprinter, a fax machine, or non-visual displays such as audio outputdevices. The display subsystem may be a cathode ray tube (CRT), aflat-panel device such as a liquid crystal display (LCD), a projectiondevice, or the like. The display subsystem may also provide a non-visualdisplay such as via audio output devices. In general, use of the term“output device” is intended to include a variety of conventional andproprietary devices and ways to output information from computer system22 to a user.

Storage subsystem 56 stores the basic programming and data constructsthat provide the functionality of the various embodiments of the presentinvention. For example, a database and modules implementing thefunctionality of the methods of the present invention, as describedherein, may be stored in storage subsystem 56. These software modulesare generally executed by processor 52. In a distributed environment,the software modules may be stored on a plurality of computer systemsand executed by processors of the plurality of computer systems. Storagesubsystem 56 typically comprises memory subsystem 58 and file storagesubsystem 60.

Memory subsystem 58 typically includes a number of memories including amain random access memory (RAM) 70 for storage of instructions and dataduring program execution and a read only memory (ROM) 72 in which fixedinstructions are stored. File storage subsystem 60 provides persistent(non-volatile) storage for program and data files, and may includetangible storage media 29 (FIG. 1) which may optionally embody wavefrontsensor data, wavefront gradients, a wavefront elevation map, a treatmentmap, and/or an ablation table. File storage subsystem 60 may include ahard disk drive, a floppy disk drive along with associated removablemedia, a Compact Digital Read Only Memory (CD-ROM) drive, an opticaldrive, DVD, CD-R, CD-RW, solid-state removable memory, and/or otherremovable media cartridges or disks. One or more of the drives may belocated at remote locations on other connected computers at other sitescoupled to computer system 22. The modules implementing thefunctionality of the present invention may be stored by file storagesubsystem 60.

Bus subsystem 54 provides a mechanism for letting the various componentsand subsystems of computer system 22 communicate with each other asintended. The various subsystems and components of computer system 22need not be at the same physical location but may be distributed atvarious locations within a distributed network. Although bus subsystem54 is shown schematically as a single bus, alternate embodiments of thebus subsystem may utilize multiple busses.

Computer system 22 itself can be of varying types including a personalcomputer, a portable computer, a workstation, a computer terminal, anetwork computer, a control system in a wavefront measurement system orlaser surgical system, a mainframe, or any other data processing system.Due to the ever-changing nature of computers and networks, thedescription of computer system 22 depicted in FIG. 2 is intended only asa specific example for purposes of illustrating one embodiment of thepresent invention. Many other configurations of computer system 22 arepossible having more or less components than the computer systemdepicted in FIG. 2.

Referring now to FIG. 3, one embodiment of a wavefront measurementsystem 30 is schematically illustrated in simplified form. In verygeneral terms, wavefront measurement system 30 is configured to senselocal slopes of a gradient map exiting the patient's eye. Devices basedon the Hartmann-Shack principle generally include a lenslet array tosample the gradient map uniformly over an aperture, which is typicallythe exit pupil of the eye. Thereafter, the local slopes of the gradientmap are analyzed so as to reconstruct the wavefront surface or map.

More specifically, one wavefront measurement system 30 includes an imagesource 32, such as a laser, which projects a source image throughoptical tissues 34 of eye E so as to form an image 44 upon a surface ofretina R. The image from retina R is transmitted by the optical systemof the eye (e.g., optical tissues 34) and imaged onto a wavefront sensor36 by system optics 37. The wavefront sensor 36 communicates signals toa computer system 22′ for measurement of the optical errors in theoptical tissues 34 and/or determination of an optical tissue ablationtreatment program. Computer 22′ may include the same or similar hardwareas the computer system 22 illustrated in FIGS. 1 and 2. Computer system22′ may be in communication with computer system 22 that directs thelaser surgery system 10, or some or all of the components of computersystem 22, 22′ of the wavefront measurement system 30 and laser surgerysystem 10 may be combined or separate. If desired, data from wavefrontsensor 36 may be transmitted to a laser computer system 22 via tangiblemedia 29, via an I/O port, via an networking connection 66 such as anintranet or the Internet, or the like.

Wavefront sensor 36 generally comprises a lenslet array 38 and an imagesensor 40. As the image from retina R is transmitted through opticaltissues 34 and imaged onto a surface of image sensor 40 and an image ofthe eye pupil P is similarly imaged onto a surface of lenslet array 38,the lenslet array separates the transmitted image into an array ofbeamlets 42, and (in combination with other optical components of thesystem) images the separated beamlets on the surface of sensor 40.Sensor 40 typically comprises a charged couple device or “CCD,” andsenses the characteristics of these individual beamlets, which can beused to determine the characteristics of an associated region of opticaltissues 34. In particular, where image 44 comprises a point or smallspot of light, a location of the transmitted spot as imaged by a beamletcan directly indicate a local gradient of the associated region ofoptical tissue.

Eye E generally defines an anterior orientation ANT and a posteriororientation POS. Image source 32 generally projects an image in aposterior orientation through optical tissues 34 onto retina R asindicated in FIG. 3. Optical tissues 34 again transmit image 44 from theretina anteriorly toward wavefront sensor 36. Image 44 actually formedon retina R may be distorted by any imperfections in the eye's opticalsystem when the image source is originally transmitted by opticaltissues 34. Optionally, image source projection optics 46 may beconfigured or adapted to decrease any distortion of image 44.

In some embodiments, image source optics 46 may decrease lower orderoptical errors by compensating for spherical and/or cylindrical errorsof optical tissues 34. Higher order optical errors of the opticaltissues may also be compensated through the use of an adaptive opticelement, such as a deformable mirror (described below). Use of an imagesource 32 selected to define a point or small spot at image 44 uponretina R may facilitate the analysis of the data provided by wavefrontsensor 36. Distortion of image 44 may be limited by transmitting asource image through a central region 48 of optical tissues 34 which issmaller than a pupil 50, as the central portion of the pupil may be lessprone to optical errors than the peripheral portion. Regardless of theparticular image source structure, it will be generally be beneficial tohave a well-defined and accurately formed image 44 on retina R.

In one embodiment, the wavefront data may be stored in a computerreadable medium 29 or a memory of the wavefront sensor system 30 in twoseparate arrays containing the x and y wavefront gradient valuesobtained from image spot analysis of the Hartmann-Shack sensor images,plus the x and y pupil center offsets from the nominal center of theHartmann-Shack lenslet array, as measured by the pupil camera 51 (FIG.3) image. Such information contains all the available information on thewavefront error of the eye and is sufficient to reconstruct thewavefront or any portion of it. In such embodiments, there is no need toreprocess the Hartmann-Shack image more than once, and the data spacerequired to store the gradient array is not large. For example, toaccommodate an image of a pupil with an 8 mm diameter, an array of a20×20 size (i.e., 400 elements) is often sufficient. As can beappreciated, in other embodiments, the wavefront data may be stored in amemory of the wavefront sensor system in a single array or multiplearrays.

While the methods of the present invention will generally be describedwith reference to sensing of an image 44, it should be understood that aseries of wavefront sensor data readings may be taken. For example, atime series of wavefront data readings may help to provide a moreaccurate overall determination of the ocular tissue aberrations. As theocular tissues can vary in shape over a brief period of time, aplurality of temporally separated wavefront sensor measurements canavoid relying on a single snapshot of the optical characteristics as thebasis for a refractive correcting procedure. Still further alternativesare also available, including taking wavefront sensor data of the eyewith the eye in differing configurations, positions, and/ororientations. For example, a patient will often help maintain alignmentof the eye with wavefront measurement system 30 by focusing on afixation target, as described in U.S. Pat. No. 6,004,313, the fulldisclosure of which is incorporated herein by reference. By varying aposition of the fixation target as described in that reference, opticalcharacteristics of the eye may be determined while the eye accommodatesor adapts to image a field of view at a varying distance and/or angles.

The location of the optical axis of the eye may be verified by referenceto the data provided from a pupil camera 52. In the exemplaryembodiment, a pupil camera 52 images pupil 50 so as to determine aposition of the pupil for registration of the wavefront sensor datarelative to the optical tissues.

An alternative embodiment of a wavefront measurement system isillustrated in FIG. 3A. The major components of the system of FIG. 3Aare similar to those of FIG. 3. Additionally, FIG. 3A includes anadaptive optical element 53 in the form of a deformable mirror. Thesource image is reflected from deformable mirror 98 during transmissionto retina R, and the deformable mirror is also along the optical pathused to form the transmitted image between retina R and imaging sensor40. Deformable mirror 98 can be controllably deformed by computer system22 to limit distortion of the image formed on the retina or ofsubsequent images formed of the images formed on the retina, and mayenhance the accuracy of the resultant wavefront data. The structure anduse of the system of FIG. 3A are more fully described in U.S. Pat. No.6,095,651, the full disclosure of which is incorporated herein byreference.

The components of an embodiment of a wavefront measurement system formeasuring the eye and ablations comprise elements of a VISX WaveScan®,available from VISX, INCORPORATED of Santa Clara, Calif. One embodimentincludes a WaveScan® with a deformable mirror as described above. Analternate embodiment of a wavefront measuring system is described inU.S. Pat. No. 6,271,915, the full disclosure of which is incorporatedherein by reference.

The use of modal reconstruction with Zernike polynomials, as well as acomparison of modal and zonal reconstructions, has been discussed indetail by W. H. Southwell, “Wave-front estimation from wave-front slopemeasurements,” J. Opt. Soc. Am. 70:998-1006 (1980). Relatedly, G. Dai,“Modal wave-front reconstruction with Zernike polynomials andKarhunen-Loeve functions,” J. Opt. Soc. Am. 13:1218-1225 (1996) providesa detailed analysis of various wavefront reconstruction errors withmodal reconstruction with Zernike polynomials. Zernike polynomials havebeen employed to model the optical surface, as proposed by Liang et al.,in “Objective Measurement of Wave Aberrations of the Human Eye with theUse of a Harman-Shack Wave-front Sensor,” J. Opt. Soc. Am.11(7):1949-1957 (1994). The entire contents of each of these referencesare hereby incorporated by reference.

The Zernike function method of surface reconstruction and its accuracyfor normal eyes have been studied extensively for regular corneal shapesin Schweigerling, J. et al., “Using corneal height maps and polynomialdecomposition to determine corneal aberrations,” Opt. Vis. Sci., Vol.74, No. 11 (1997) and Guirao, A. et al. “Corneal wave aberration fromvideokeratography: Accuracy and limitations of the procedure,” J. Opt.Soc. Am., Vol. 17, No. 6 (2000). D. R. Ishkander et al., “An AlternativePolynomial Representation of the Wavefront Error Function,” IEEETransactions on Biomedical Engineering, Vol. 49, No. 4, (2002) reportthat the 6th order Zernike polynomial reconstruction method provides aninferior fit when compared to a method of Bhatia-Wolf polynomials. Theentire contents of each of these references are hereby incorporated byreference.

Modal wavefront reconstruction typically involves expanding thewavefront into a set of basis functions. Use of Zernike polynomials as awavefront expansion basis function has been accepted in the wavefronttechnology field due to the fact that Zernike polynomials are a set ofcomplete and orthogonal functions over a circular pupil. In addition,some lower order Zernike modes, such as defocus, astigmatism, coma andspherical aberrations, represent classical aberrations. Unfortunately,there may be drawbacks to the use of Zernike polynomials. Because theZernike basis function has a rapid fluctuation near the aperture,especially for higher orders, a slight change in the Zernikecoefficients can greatly affect the wavefront surface. Further, due tothe aberration balancing between low and high order Zernike modes,truncation of Zernike series often causes inconsistent Zernikecoefficients.

In order to solve some of the above-mentioned problems with Zernikereconstruction, we looked to other basis functions. Fourier seriesappear to be an advantageous basis function set due to its robust fastFourier transform (FFT) algorithm. Also, the derivatives of Fourierseries are still a Fourier series. For un-bounded functions (i.e. withno boundary conditions), Fourier reconstruction can be used to directlyestimate the function from a set of gradient data. It may be difficult,however, to apply this technique directly to wavefront technologybecause wavefront reconstruction typically relates to a boundedfunction, or a function with a pupil aperture.

Iterative Fourier reconstruction techniques can apply to boundedfunctions with unlimited aperture functions. This is to say, theaperture of the function can be circular, annular, oval, square,rectangular, or any other shape. Such an approach is discussed inRoddier et al., “Wavefront reconstruction using iterative Fouriertransforms,” Appl. Opt. 30, 1325-1327 (1991), the entire contents ofwhich are hereby incorporated by reference. Such approaches, however,are significantly improved by accounting for missing data points due tocorneal reflection, bad CCD pixels, and the like.

I. Determining an Optical Surface Model for an Optical Tissue System ofan Eye

The present invention provides systems, software, and methods that use ahigh speed and accurate iterative Fourier transformation algorithm tomathematically determine an optical surface model for an optical tissuesystem of an eye.

A. Inputting Optical Data from the Optical Tissue System of the Eye

There are a variety of devices and methods for generating optical datafrom optical tissue systems. The category of aberroscopes oraberrometers includes classical phoropter and wavefront approaches.Topography based measuring devices and methods can also be used togenerate optical data. Wavefront devices are often used to measure bothlow order and high order aberrations of an optical tissue system.Particularly, wavefront analysis typically involves transmitting animage through the optical system of the eye, and determining a set ofsurface gradients for the optical tissue system based on the transmittedimage. The surface gradients can be used to determine the optical data.

B. Determining the Optical Surface Model by Applying an IterativeFourier Transform to the Optical Data

FIG. 4 schematically illustrates a simplified set of modules forcarrying out a method according to one embodiment of the presentinvention. The modules may be software modules on a computer readablemedium that is processed by processor 52 (FIG. 2), hardware modules, ora combination thereof. A wavefront aberration module 80 typicallyreceives data from the wavefront sensors and measures the aberrationsand other optical characteristics of the entire optical tissue systemimaged. The data from the wavefront sensors are typically generated bytransmitting an image (such as a small spot or point of light) throughthe optical tissues, as described above. Wavefront aberration module 80produces an array of optical gradients or a gradient map. The opticalgradient data from wavefront aberration module 80 may be transmitted toa Fourier transform module 82, where an optical surface or a modelthereof, or a wavefront elevation surface map, can be mathematicallyreconstructed from the optical gradient data.

It should be understood that the optical surface or model thereof neednot precisely match an actual tissue surface, as the gradient data willshow the effects of aberrations which are actually located throughoutthe ocular tissue system. Nonetheless, corrections imposed on an opticaltissue surface so as to correct the aberrations derived from thegradients should correct the optical tissue system. As used herein termssuch as “an optical tissue surface” or “an optical surface model” mayencompass a theoretical tissue surface (derived, for example, fromwavefront sensor data), an actual tissue surface, and/or a tissuesurface formed for purposes of treatment (for example, by incisingcorneal tissues so as to allow a flap of the corneal epithelium andstroma to be displaced and expose the underlying stroma during a LASIKprocedure).

Once the wavefront elevation surface map is generated by Fouriertransform module 82, the wavefront gradient map may be transmitted to alaser treatment module 84 for generation of a laser ablation treatmentto treat or ameliorate optical errors in the optical tissues.

FIG. 5 is a detailed flow chart which illustrates a data flow and methodsteps of one Fourier based method of generating a laser ablationtreatment. The illustrated method is typically carried out by a systemthat includes a processor and a memory coupled to the processor. Thememory may be configured to store a plurality of modules which have theinstructions and algorithms for carrying out the steps of the method.

As can be appreciated, the present invention should not be limited tothe order of steps, or the specific steps illustrated, and variousmodifications to the method, such as having more or less steps, may bemade without departing from the scope of the present invention. Forcomparison purposes, a series expansion method of generating a wavefrontelevation map is shown in dotted lines, and are optional steps.

A wavefront measurement system that includes a wavefront sensor (such asa Hartmann-Shack sensor) may be used to obtain one or more displacementmaps 90 (e.g., Hartmann-Shack displacement maps) of the optical tissuesof the eye. The displacement map may be obtained by transmitting animage through the optical tissues of the eye and sensing the exitingwavefront surface.

From the displacement map 90, it is possible to calculate a surfacegradient or gradient map 92 (e.g., Hartmann-Shack gradient map) acrossthe optical tissues of the eye. Gradient map 92 may comprise an array ofthe localized gradients as calculated from each location for eachlenslet of the Hartmann-Shack sensor.

A Fourier transform may be applied to the gradient map to mathematicallyreconstruct the optical tissues or to determine an optical surfacemodel. The Fourier transform will typically output the reconstructedoptical tissue or the optical surface model in the form of a wavefrontelevation map. For the purposes of the instant invention, the termFourier transform also encompasses iterative Fourier transforms.

It has been found that a Fourier transform reconstruction method, suchas a fast Fourier transformation (FFT), is many times faster thancurrently used 6th order Zernike or polynomial reconstruction methodsand yields a more accurate reconstruction of the actual wavefront.Advantageously, the Fourier reconstruction limits the specialfrequencies used in reconstruction to the Nyquist limit for the datadensity available and gives better resolution without aliasing. If it isdesired, for some a priori reason, to limit the spatial frequenciesused, this can be done by truncating the transforms of the gradient inFourier transformation space midway through the calculation. If it isdesired to sample a small portion of the available wavefront or decenterit, this may be done with a simple mask operation on the gradient databefore the Fourier transformation operation. Unlike Zernikereconstruction methods in which the pupil size and centralization of thepupil is required, such concerns do not effect the fast Fouriertransformation.

Moreover, since the wavefront sensors measure x- and y-components of thegradient map on a regularly spaced grid, the data is band-limited andthe data contains no spatial frequencies larger than the Nyquist ratethat corresponds to the spacing of the lenslets in the instrument(typically, the lenslets will be spaced no more than about 0.8 mm andabout 0.1 mm, and typically about 0.4 mm). Because the data is on aregularly spaced Cartesian grid, non-radial reconstruction methods, suchas a Fourier transform, are well suited for the band-limited data.

In contrast to the Fourier transform, when a series expansion techniqueis used to generate a wavefront elevation map 100 from the gradient map92, the gradient map 92 and selected expansion series 96 are used toderive appropriate expansion series coefficients 98. A particularlybeneficial form of a mathematical series expansion for modeling thetissue surface are Zernike polynomials. Typical Zernike polynomial setsincluding terms 0 through 6th order or 0 through 10th order are used.The coefficients an for each Zernike polynomial Z_(n) may, for example,be determined using a standard least squares fit technique. The numberof Zernike polynomial coefficients a_(n) may be limited (for example, toabout 28 coefficients).

While generally considered convenient for modeling of the opticalsurface so as to generate an elevation map, Zernike polynomials (andperhaps all series expansions) can introduce errors. Nonetheless,combining the Zernike polynomials with their coefficients and summingthe Zernike coefficients 99 allows a wavefront elevation map 100 to becalculated, and in some cases, may very accurately reconstruct awavefront elevation map 100.

It has been found that in some instances, especially where the error inthe optical tissues of the eye is spherical, the Zernike reconstructionmay be more accurate than the Fourier transform reconstruction. Thus, insome embodiments, the modules of the present invention may include botha Fourier transform module 94 and Zernike modules 96, 98, 99. In suchembodiments, the reconstructed surfaces obtained by the two modules maybe compared by a comparison module (not shown) to determine which of thetwo modules provides a more accurate wavefront elevation map. The moreaccurate wavefront elevation map may then be used by 100, 102 tocalculate the treatment map and ablation table, respectively.

In one embodiment, the wavefront elevation map module 100 may calculatethe wavefront elevation maps from each of the modules and a gradientfield may be calculated from each of the wavefront elevation maps. Inone configuration, the comparison module may apply a merit function todetermine the difference between each of the gradient maps and anoriginally measured gradient map. One example of a merit function is theroot mean square gradient error, RMS_(grad), found from the followingequation:${RMS}_{grad} = \sqrt{\frac{\sum\limits_{alldatapoints}\begin{matrix}\left\{ {\left( {\frac{\partial{W\left( {x,y} \right)}}{\partial x} - {{Dx}\left( {x,y} \right)}^{2}} \right) +} \right. \\\left. \left( {\frac{\partial{W\left( {x,y} \right)}}{\partial y} - {{Dy}\left( {x,y} \right)}^{2}} \right) \right\}\end{matrix}}{N}}$

-   -   where:    -   N is the number of locations sampled    -   (x,y) is the sample location    -   δW(x,y)/δx is the x component of the reconstructed wavefront        gradient    -   δW(x,y)/δy is the y component of the reconstructed wavefront        gradient    -   Dx(x,y) is the x component of the gradient data    -   Dy(x,y) is the y component of the gradient data

If the gradient map from the Zernike reconstruction is more accurate,the Zernike reconstruction is used. If the Fourier reconstruction ismore accurate, the Fourier reconstruction is used.

After the wavefront elevation map is calculated, treatment map 102 maythereafter be calculated from the wavefront elevation map 100 so as toremove the regular (spherical and/or cylindrical) and irregular errorsof the optical tissues. By combining the treatment map 102 with a laserablation pulse characteristics 104 of a particular laser system, anablation table 106 of ablation pulse locations, sizes, shapes, and/ornumbers can be developed.

A laser treatment ablation table 106 may include horizontal and verticalposition of the laser beam on the eye for each laser beam pulse in aseries of pulses. The diameter of the beam may be varied during thetreatment from about 0.65 mm to 6.5 mm. The treatment ablation table 106typically includes between several hundred pulses to five thousand ormore pulses, and the number of laser beam pulses varies with the amountof material removed and laser beam diameters employed by the lasertreatment table. Ablation table 106 may optionally be optimized bysorting of the individual pulses so as to avoid localized heating,minimize irregular ablations if the treatment program is interrupted,and the like. The eye can thereafter be ablated according to thetreatment table 106 by laser ablation.

In one embodiment, laser ablation table 106 may adjust laser beam 14 toproduce the desired sculpting using a variety of alternative mechanisms.The laser beam 14 may be selectively limited using one or more variableapertures. An exemplary variable aperture system having a variable irisand a variable width slit is described in U.S. Pat. No. 5,713,892, thefull disclosure of which is incorporated herein by reference. The laserbeam may also be tailored by varying the size and offset of the laserspot from an axis of the eye, as described in U.S. Pat. Nos. 5,683,379and 6,203,539, and as also described in U.S. application Ser. No.09/274,999 filed Mar. 22, 1999, the full disclosures of which areincorporated herein by reference.

Still further alternatives are possible, including scanning of the laserbeam over a surface of the eye and controlling the number of pulsesand/or dwell time at each location, as described, for example, by U.S.Pat. No. 4,665,913 (the full disclosure of which is incorporated hereinby reference); using masks in the optical path of laser beam 14 whichablate to vary the profile of the beam incident on the cornea, asdescribed in U.S. patent application Ser. No. 08/468,898, filed Jun. 6,1995 (the full disclosure of which is incorporated herein by reference);hybrid profile-scanning systems in which a variable size beam (typicallycontrolled by a variable width slit and/or variable diameter irisdiaphragm) is scanned across the cornea; or the like.

One exemplary method and system for preparing such an ablation table isdescribed in U.S. Pat. No. 6,673,062, the full disclosure of which isincorporated herein by reference.

The mathematical development for the surface reconstruction from surfacegradient data using a Fourier transform algorithm according to oneembodiment of the present invention will now be described. Suchmathematical algorithms will typically be incorporated by Fouriertransform module 82 (FIG. 4), Fourier transform step 94 (FIG. 5), orother comparable software or hardware modules to reconstruct thewavefront surface. As can be appreciated, the Fourier transformalgorithm described below is merely an example, and the presentinvention should not be limited to this specific implementation.

First, let there be a surface that may be represented by the functions(x,y) and let there be data giving the gradients of this surface,$\frac{\partial{s\left( {x,y} \right)}}{\partial x}\quad{and}\quad{\frac{\partial{s\left( {x,y} \right)}}{\partial y}.}$The goal is to find the surface s(x,y) from the gradient data.

Let the surface be locally integratable over all space so that it may berepresented by a Fourier transform. The Fourier transform of the surfaceis then given by $\begin{matrix}{{F\left( {s\left( {x,y} \right)} \right)} = {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{s\left( {x,y} \right)}{\mathbb{e}}^{- {{\mathbb{i}}{({{ux} + {vy}})}}}{\mathbb{d}x}{\mathbb{d}y}}}}} = {S\left( {u,v} \right)}}} & (1)\end{matrix}$

The surface may then be reconstructed from the transform coefficients,S(u,v), using $\begin{matrix}{{s\left( {x,y} \right)} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{S\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}} & (2)\end{matrix}$

Equation (2) may now be used to give a representation of the x componentof the gradient, $\frac{\partial{s\left( {x,y} \right)}}{\partial x}$in terms of the Fourier coefficients for the surface:$\frac{\partial{s\left( {x,y} \right)}}{\partial x} = \frac{\partial\left( {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{S\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} \right)}{\partial x}$

Differentiation under the integral then gives: $\begin{matrix}{\frac{\partial{s\left( {x,y} \right)}}{\partial x} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{uS}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}} & (3)\end{matrix}$

A similar equation to (3) gives a representation of the y component ofthe gradient in terms of the Fourier coefficients:$\frac{\partial{s\left( {x,y} \right)}}{\partial y} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{vS}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}$

The x component of the gradient is a function of x and y so it may alsobe represented by the coefficients resulting from a Fouriertransformation. Let the${{dx}\left( {x,y} \right)} = \frac{\partial{s\left( {x,y} \right)}}{\partial x}$so that, following the logic that led to (2) $\begin{matrix}{{{{dx}\left( {x,y} \right)} = {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{{Dx}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} = \frac{\partial{s\left( {x,y} \right)}}{\partial x}}}{where}} & (5) \\{{F\left( {{dx}\left( {x,y} \right)} \right)} = {{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{{dx}\left( {x,y} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}} = {{Dx}\left( {u,v} \right)}}} & (6)\end{matrix}$

Equation (3) must equal (5) and by inspecting similar terms it may beseen that in general this can only be true if

-   -   Dx(u,v)=iuS(u,v)        or $\begin{matrix}        {{S\left( {u,v} \right)} = \frac{- {{iDx}\left( {u,v} \right)}}{u}} & (7)        \end{matrix}$

A similar development for the y gradient component using (4) leads to$\begin{matrix}{{S\left( {u,v} \right)} = \frac{- {{iDy}\left( {u,v} \right)}}{v}} & (8)\end{matrix}$

Note that (7) and (8) indicate that Dx (u,v) and Dy(u,v) arefunctionally dependent with the relationshipvDx(u,v)=uDy(u,v)

The surface may now be reconstructed from the gradient data by firstperforming a discrete Fourier decomposition of the two gradient fields,dx and dy to generate the discrete Fourier gradient coefficients Dx(u,v)and Dy(u,v). From these components (7) and (8) are used to find theFourier coefficients of the surface S(u,v). These in turn are used withan inverse discrete Fourier transform to reconstruct the surface s(x,y).

The above treatment makes a non-symmetrical use of the discrete Fouriergradient coefficients in that one or the other is used to find theFourier coefficients of the surface. The method makes use of theLaplacian, a polynomial, second order differential operator, given by$L \equiv {\frac{\partial^{2}}{\partial^{2}x} + {\frac{\partial^{2}}{\partial^{2}y}\quad{or}\quad L}} \equiv {{\frac{\partial}{\partial x}\quad\left( \frac{\partial}{\partial x} \right)} + {\frac{\partial}{\partial y}\quad\left( \frac{\partial}{\partial y} \right)}}$

So when the Laplacian acts on the surface function, s(x,y), the resultis${{Ls}\left( {x,y} \right)} = {\frac{\partial^{2}{s\left( {x,y} \right)}}{\partial^{2}x} + {\frac{\partial^{2}{s\left( {x,y} \right)}}{\partial^{2}y}\quad{or}}}$${{Ls}\left( {x,y} \right)} = {{\frac{\partial}{\partial x}\left( \frac{\partial{s\left( {x,y} \right)}}{\partial x} \right)} + {\frac{\partial}{\partial y}\left( \frac{\partial{s\left( {x,y} \right)}}{\partial y} \right)}}$

Using the second form given above and substituting (3) for$\frac{\partial{s\left( {x,y} \right)}}{\partial x}$in the first integral of the sum and (4) for$\frac{\partial{s\left( {x,y} \right)}}{\partial y}$in the second term, the Laplacian of the surface is found to be$\begin{matrix}\begin{matrix}{{{Ls}\left( {x,y} \right)} = {{\frac{\partial}{\partial x}\left( {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{uS}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} \right)} +}} \\{\frac{\partial}{\partial y}\left( {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{vS}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} \right)}\end{matrix} & \quad \\\begin{matrix}{{{Ls}\left( {x,y} \right)} = {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{- u^{2}}{S\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} +}} \\{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{- v^{2}}{S\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}\end{matrix} & \quad \\{{{Ls}\left( {x,y} \right)} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{- \left( {u^{2} + v^{2}} \right)}{S\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}} & (9)\end{matrix}$

Equation (9) shows that the Fourier coefficients of the Laplacian of atwo dimensional function are equal to −(u²+v²) times the Fouriercoefficients of the function itself so that${S\left( {u,v} \right)} = \frac{- {F\left( {{Ls}\left( {x,y} \right)} \right)}}{\left( {u^{2} + v^{2}} \right)}$

Now let the Laplacian be expressed in terms of dx(x,y) and dy(x,y) asdefined above so that${{Ls}\left( {x,y} \right)} = {{\frac{\partial}{\partial x}\left( {{dx}\left( {x,y} \right)} \right)} + {\frac{\partial}{\partial y}\left( {{dy}\left( {x,y} \right)} \right)}}$and through the use of (5) and the similar expression for dy(x,y)$\begin{matrix}\begin{matrix}{{{Ls}\left( {x,y} \right)} = {{\frac{\partial}{\partial x}\left( {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{{Dx}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} \right)} +}} \\{\frac{\partial}{\partial y}\left( {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{{Dy}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} \right)}\end{matrix} & \quad \\\begin{matrix}{{{Ls}\left( {x,y} \right)} = {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{uDx}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}} +}} \\{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{i}}\quad{{vDy}\left( {u,v} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}\end{matrix} & \quad \\{{{Ls}\left( {x,y} \right)} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{{\mathbb{i}}\left( {{{uDx}\left( {u,v} \right)} + {{vDy}\left( {u,v} \right)}} \right)}{\mathbb{e}}^{{\mathbb{i}}{({{ux} + {vy}})}}{\mathbb{d}u}{\mathbb{d}v}}}}}} & (10)\end{matrix}$(9) and (10) must be equal and comparing them, it is seen that thisrequires that: $\begin{matrix}{{{{- \left( {u^{2} + v^{2}} \right)}{S\left( {u,v} \right)}} = {{i\left( {{{uDx}\left( {u,v} \right)} + {{vD}\left( {u,v} \right)}} \right)}\quad{or}}}{{S\left( {u,v} \right)} = \frac{- {i\left( {{{uDx}\left( {u,v} \right)} + {{vDy}\left( {u,v} \right)}} \right)}}{u^{2} + v^{2}}}} & (11)\end{matrix}$

As before, Dx(u,v) and Dy(u,v) are found by taking the Fouriertransforms of the measured gradient field components. They are then usedin (11) to find the Fourier coefficients of the surface itself, which inturn is reconstructed from them. This method has the effect of using allavailable information in the reconstruction, whereas the Zernikepolynomial method fails to use all of the available information.

It should be noted, s(x,y) may be expressed as the sum of a new functions(x,y)′ and a tilted plane surface passing through the origin. This sumis given by the equations(x,y)=s(x,y)′+ax+by

Then the first partial derivatives of f(x,y) with respect to x and y aregiven by$\frac{\partial{s\left( {x,y} \right)}}{\partial x} = {\frac{\partial{s\left( {x,y} \right)}^{\prime}}{\partial x} + a}$$\frac{\partial{s\left( {x,y} \right)}}{\partial y} = {\frac{\partial{s\left( {x,y} \right)}^{\prime}}{\partial y} + b}$

Now following the same procedure that lead to (6), the Fourier transformof these partial derivatives, Dx(u,v) and Dy(u,v), are found to be$\begin{matrix}\left. {\begin{matrix}{{{Dx}\left( {u,v} \right)} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\frac{\partial{s\left( {x,y} \right)}}{\partial x}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}}} \\{= {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\frac{\partial{s\left( {x,y} \right)}^{\prime}}{\partial x}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}} +}} \\{\frac{a}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}}\end{matrix}{{{Dx}\left( {u,v} \right)} = {{{Dx}\left( {u,v} \right)}^{\prime} + {a\quad{\delta\left( {u,v} \right)}}}}} \right) & (12) \\\left. {\begin{matrix}{{{Dy}\left( {u,v} \right)} = {\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\frac{\partial{s\left( {x,y} \right)}}{\partial y}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}}} \\{= {{\frac{1}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\frac{\partial{s\left( {x,y} \right)}^{\prime}}{\partial y}{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}} +}} \\{\frac{b}{2\quad\pi}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{e}}^{{- {\mathbb{i}}}\quad{({{wx} + {vy}})}}{\mathbb{d}x}{\mathbb{d}y}}}}}\end{matrix}{{{Dy}\left( {u,v} \right)} = {{{Dy}\left( {u,v} \right)}^{\prime} + {b\quad{\delta\left( {u,v} \right)}}}}} \right) & (13)\end{matrix}$

In (12) and (13), δ(u,v), is the Dirac delta function that takes thevalue 1 if u=v=0 and takes the value 0 otherwise. Using (12) and (13) in(11), the expression of the Fourier transform of the surface may bewritten as${S\left( {u,v} \right)} = \frac{- {i\left( {{{uDx}\left( {u,v} \right)}^{\prime} + {{vDy}\left( {u,v} \right)}^{\prime} + {\left( {{ua} + {vb}} \right){\delta\left( {u,v} \right)}}} \right)}}{u^{2} + v^{2}}$

But it will be realized that the term in the above equation can have noeffect whatsoever on the value of S(u,v) because if u and v are not bothzero, the delta function is zero so the term vanishes. However in theonly other case, u and v are both zero and this also causes the term tovanish. This means that the surface reconstructed will not be unique butwill be a member of a family of surfaces, each member having a differenttilted plane (or linear) part. Therefore to reconstruct the uniquesurface represented by a given gradient field, the correct tilt must berestored. The tilt correction can be done in several ways.

Since the components of the gradient of the tilt plane, a and b, are thesame at every point on the surface, if the correct values can be foundat one point, they can be applied everywhere at once by adding to theinitially reconstructed surface, s(x,y)′, the surface (ax+by). This maybe done by finding the gradient of s(x,y)′ at one point and subtractingthe components from the given components. The differences are theconstant gradient values a and b. However when using real data there isalways noise and so it is best to use an average to find a and b. Oneuseful way to do this is to find the average gradient components of thegiven gradient field and use these averages as a and b.

δs/δx=a δs/δy=b

The reconstructed surface is then given bys(x,y)=s(x,y)′+

δs/δx

x+

δs/δy

y  (14)where s(x,y)′ is found using the Fourier reconstruction method developedabove.

Attention is now turned to implementation of this method using discretefast Fourier transform methods. The algorithm employed for the discretefast Fourier transform in two dimensions is given by the equation$\begin{matrix}{{F\left( {k,l} \right)} = {\sum\limits_{m = 1}^{M}{\sum\limits_{n = 1}^{N}{{f\left( {n,m} \right)}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{{({k - 1})}{({n - 1})}}{N} + \frac{{({l - 1})}{({m - 1})}}{M}})}}}}}}} & (12)\end{matrix}$with the inverse transform given by $\begin{matrix}{{f\left( {n,m} \right)} = {\frac{1}{MN}{\sum\limits_{k = 1}^{M}{\sum\limits_{l = 1}^{N}{{F\left( {k,l} \right)}{\mathbb{e}}^{{\mathbb{i}}\quad 2\quad{\pi{({\frac{{({k - 1})}{({n - 1})}}{N} + \frac{{({l - 1})}{({m - 1})}}{M}})}}}}}}}} & (13)\end{matrix}$

When these equations are implemented N and M are usually chosen so to beequal. For speed of computation, they are usually taken to be powers of2. (12) and (13) assume that the function is sampled at locationsseparated by intervals dx and dy. For reasons of algorithmicsimplification, as shown below, dx and dy are usually set equal.

In equation (12) let n be the index of the x data in array f(n,m) andlet k be the index of the variable u in the transform array, F(k,l).

Let us begin by supposing that in the discrete case the x values arespaced by a distance dx, so when equation (12) is used, each time n isincremented, x is changed by an amount dx. So by choosing the coordinatesystem properly we could represent the position of the pupil data x by:x=(n−1)dxso that:(n−1)=x/dx

Likewise, (m−1) may be set equal to y/dy. Using these relationships,(12) may be written as:${F\left( {k,l} \right)} = {\sum\limits_{m = 1}^{M}{\sum\limits_{n = 1}^{N}{{f\left( {n,m} \right)}{\mathbb{e}}^{{- {\mathbb{i}}}\quad 2\quad{\pi{({\frac{{({k - 1})}x}{Ndx} + \frac{{({l - 1})}y}{Mdy}})}}}}}}$

Comparing the exponential terms in this discrete equation with those inits integral form (1), it is seen that$u = {{\frac{2\quad{\pi\left( {k - 1} \right)}}{Ndx}\quad{and}\quad v} = \frac{2\quad{\pi\left( {l - 1} \right)}}{Mdy}}$

In these equations notice that Ndx is the x width of the sampled areaand Mdy is the y width of the sampled area. LettingNdx=X, the total x widthMdy=Y, the total y width  (14)

The above equations become $\begin{matrix}{{u(k)} = {{\frac{2\quad{\pi\left( {k - 1} \right)}}{X}\quad{and}\quad{v(l)}} = \frac{2\quad{\pi\left( {l - 1} \right)}}{Y}}} & (15)\end{matrix}$

Equations (15) allow the Fourier coefficients, Dx(k,l) and Dy(k,l),found from the discrete fast Fourier transform of the gradientcomponents, dx(n,m) and dy(n,m), to be converted into the discreteFourier coefficients of the surface, S(k,l) as follows.${S\left( {k,l} \right)} = \frac{{{- {\mathbb{i}}}\frac{2\quad{\pi\left( {k - 1} \right)}}{X}{{Dx}\left( {k,l} \right)}} - {{\mathbb{i}}\frac{2\quad{\pi\left( {l - 1} \right)}}{Y}{{Dy}\left( {k,l} \right)}}}{\left( \frac{2\quad{\pi\left( {k - 1} \right)}}{X} \right)^{2} + \left( \frac{2\quad{\pi\left( {l - 1} \right)}}{Y} \right)^{2}}$

This equation is simplified considerably if the above mentioned N ischosen equal to M and dx is chosen to dy, so that X=Y. It then becomes$\begin{matrix}{{S\left( {k,1} \right)} = {\left( \frac{{- i}\quad X}{2\quad\pi} \right)\frac{{\left( {k - 1} \right){{Dx}\left( {k,1} \right)}} + {\left( {1 - 1} \right){{Dy}\left( {k,1} \right)}}}{\left( {k - 1} \right)^{2} + \left( {1 - 1} \right)^{2}}}} & (16)\end{matrix}$

Let us now consider the values u(k) and v(l) take in (13) as k variesfrom 1 to N and l varies from 1 to M. When k=l=1, u=v=0 and theexponential takes the value 1. As u and v are incremented by 1 so thatk=l=2, u and v are incremented by unit increments, du and dv, so${u(2)} = {\frac{2\quad\pi}{X} = {{{du}\quad{and}\quad{v(2)}} = {\frac{2\quad\pi}{Y} = {dv}}}}$

Each increment of k or l increments u or v by amounts du and dvrespectively so that for any value of k or l, u(k) and v(l) may bewritten asu(k)=(k−1)du, v(l)=(l−1)dv

This process may be continued until k=N and l=M at which time${u(N)} = {\frac{2\quad{\pi\left( {N - 1} \right)}}{X} = {{\frac{2\quad{\pi(N)}}{X} - \frac{2\quad\pi}{X}} = {\frac{2\quad{\pi(N)}}{X} - {du}}}}$${v(M)} = {\frac{2\quad{\pi\left( {M - 1} \right)}}{X} = {{\frac{2\quad{\pi(M)}}{X} - \frac{2\quad\pi}{X}} = {\frac{2\quad{\pi(M)}}{X} - {dv}}}}$

But now consider the value the exponential in (13) takes when theseconditions hold. In the following, the exponential is expressed as aproduct $\begin{matrix}{{{\mathbb{e}}^{\frac{{\mathbb{i}}\quad 2\quad{\pi{({N - 1})}}{({n - 1})}}{N}}{\mathbb{e}}^{\frac{{\mathbb{i}}\quad 2\quad{\pi{({M - 1})}}{({m - 1})}}{M}}} = {{\mathbb{e}}^{({{{\mathbb{i}}\quad 2{\pi{({n - {1c}})}}} - \frac{{\mathbb{i}}\quad 2\quad{\pi{({n - 1})}}}{N}})}{\mathbb{e}}^{({{{\mathbb{i}}\quad 2\quad{\pi{({m - {1c}})}}} - \frac{{\mathbb{i}}\quad 2\quad{\pi{({m - 1})}}}{M}})}}} \\{= {{\mathbb{e}}^{\frac{{- {\mathbb{i}}}\quad 2\quad{\pi{({n - 1})}}}{N}}{\mathbb{e}}^{\frac{{- {\mathbb{i}}}\quad 2\quad{\pi{({m - 1})}}}{M}}}}\end{matrix}$

Using the same logic as was used to obtain equations (15), the values ofu(N) and v(M) areu(N)=−du and v(M)=−dv

Using this fact, the following correlation may now be made betweenu(k+1) and u(N−k) and between v(l+1) and v(M−1) for k>1 and l>1u(k)=−u(N−k+2) v(l)=−v(M−l+2)

In light of equations (15) $\begin{matrix}{{u\left( {N - k + 2} \right)} = {{\frac{{- 2}\quad{\pi\left( {k - 1} \right)}}{X}\quad{and}\quad{v\left( {M - 1 + 2} \right)}} = \frac{{- 2}\quad{\pi\left( {1 - 1} \right)}}{Y}}} & (17)\end{matrix}$

To implement (15), first note that Dx(k,l) and Dy(k,l) are formed asmatrix arrays and so it is best to form the coefficients (k−1) and (l−1)as matrix arrays so that matrix multiplication method may be employed toform S(k,l) as a matrix array.

Assuming the Dx and Dy are square arrays of N×N elements, let the (k−1)array, K(k,l) be formed as a N×N array whose rows are all the sameconsisting of integers starting at 0 (k=1) and progressing with integerinterval of 1 to the value ceil(N/2)−1. The “ceil” operator rounds thevalue N/2 to the next higher integer and is used to account for caseswhere N is an odd number. Then, in light of the relationships given in(17), the value of the next row element is given the value −floor(N/2).The “floor” operator rounds the value N/2 to the next lower integer,used again for cases where N is odd. The row element following the onewith value −floor(N/2) is incremented by 1 and this proceeds until thelast element is reached (k=N) and takes the value −1. In this way, whenmatrix |Dx| is multiplied term by term times matrix |K|, each term of|Dx| with the same value of k is multiplied by the correct integer andhence by the correct u value.

Likewise, let matrix |L(k,l)| be formed as an N×N array whose columnsare all the same consisting of integers starting at 0 (l=1) andprogressing with integer interval of 1 to the value ceil(N/2)−1. Then,in light of the relationships given in (17), the value of the nextcolumn element is given the value −floor(N/2). The column elementfollowing the one with value −floor(N/2) is incremented by 1 and thisproceeds until the last element is reached (l=N) and takes the value −1.In this way, when matrix |Dy| is multiplied term by term times matrix|L|, each term of |Dy| with the same value of 1 is multiplied by thecorrect integer and hence by the correct v value.

The denominator of (15) by creating a matrix |D| from the sum ofmatrices formed by multiplying, term-by-term, |K| times itself and |L|times itself. The (1,1) element of |D| is always zero and to avoiddivide by zero problems, it is set equal to 1 after |D| is initiallyformed. Since the (1,1) elements of |K| and |L| are also zero at thistime, this has the effect of setting the (1,1) element of |S| equal tozero. This is turn means that the average elevation of the reconstructedsurface is zero as may be appreciated by considering that the value of(12) when k=l=1 is the sum of all values of the function f(x,y). If thissum is zero, the average value of the function is zero.

Let the term-by-term multiplication of two matrices |A| and |B| besymbolized by |A|.*|B| and the term-by-term division of |A| by |B| by|A|./|B|. Then in matrix form, (16) may be written as: $\begin{matrix}{{S} = {\left( \frac{{- i}\quad X}{2\quad\pi} \right){{\left( {{{{K}.}*{{Dx}}} + {{{L}.}*{{Dy}\left( {k,1} \right)}}} \right).}/\left( {{{{K}.}*{K}} + {{{L}.}*{L}}} \right)}}} & (18)\end{matrix}$

The common factor $\left( \frac{{- i}\quad X}{2\quad\pi} \right)$is neither a function of position nor “frequency” (the variables of theFourier transform space). It is therefore a global scaling factor.

As a practical matter when coding (18), it is simpler to form K and L ifthe transform matrices Dx and Dy are first “shifted” using the standarddiscrete Fourier transform quadrant shift technique that places u=v=0element at location (floor(N/2)+1,floor(N/2)+1). The rows of K and thecolumns of L may then be formed fromrow=[1,2,3, . . . N−2,N−1,N]−(floor(N/2)+1)column=row^(T)

After the matrix |S| found with (18) using the shifted matrices, |S| isthen inverse shifted before the values of s(x,y) are found using theinverse discrete inverse Fourier transform (13).

The final step is to find the mean values of the gradient fields dx(n,m)and dy(n,m). These mean values are multiplied by the respective x and yvalues for each surface point evaluated and added to the value of s(x,y)found in the step above to give the fully reconstructed surface.

Experimental Results

A detailed description of some test methods to compare the surfacereconstructions of the expansion series (e.g., Zernike polynomial)reconstruction methods, direct integration reconstruction methods, andFourier transform reconstruction methods will now be described.

While not described in detail herein, it should be appreciated that thepresent invention also encompasses the use of direct integrationalgorithms and modules for reconstructing the wavefront elevation map.The use of Fourier transform modules, direct integration modules, andZernike modules are not contradictory or mutually exclusive, and may becombined, if desired. For example, the modules of FIG. 5 may alsoinclude direct integration modules in addition to or alternative to themodules illustrated. A more complete description of the directintegration modules and methods are described in co-pending U.S. patentapplication Ser. No. 10/006,992, filed Dec. 6, 2001 and PCT ApplicationNo. PCT/US01/46573, filed Nov. 6, 2001, both entitled “DirectWavefront-Based Corneal Ablation Treatment Program,” the completedisclosures of which are incorporated herein by reference.

To compare the various methods, a surface was ablated onto plastic, andthe various reconstruction methods were compared to a direct surfacemeasurement to determine the accuracy of the methods. Three differenttest surfaces were created for the tests, as follows:

-   -   (1) +2 ablation on a 6 mm pupil, wherein the ablation center was        offset by approximately 1 mm with respect to the pupil center;    -   (2) Presbyopia Shape I which has a 2.5 mm diameter “bump,” 1.5        μm high, decentered by 1.0 mm.    -   (3) Presbyopia Shape II which has a 2.0 mm diameter “bump,” 1.0        μm high, decentered by 0.5 mm.

The ablated surfaces were imaged by a wavefront sensor system 30 (seeFIGS. 3 and 3A), and the Hartmann-Shack spot diagrams were processed toobtain the wavefront gradients. The ablated surfaces were also scannedby a surface mapping interferometer Micro XCAM, manufactured by PhaseShift Technologies, so as to generate a high precision surface elevationmap. The elevation map directly measured by the Micro XCAM was comparedto the elevation map reconstructed by each of the different algorithms.The algorithm with the lowest root mean square (RMS) error wasconsidered to be the most effective in reconstructing the surface.

In both the direct measurement and mathematical reconstruction, theremay be a systematic “tilt” that needs correction. For the directmeasurement, the tilt in the surface (that was introduced by a tilt in asample stage holding the sample) was removed from the data bysubtracting a plane that would fit to the surface.

For the mathematical reconstructions, the angular and spatial positionsof the surface relative to the lenslet array in the wavefrontmeasurement system introduced a tilt and offset of center in thereconstruction surfaces. Correcting the “off-center” alignment was doneby identifying dominant features, such as a top of a crest, andtranslating the entire surface data to match the position of thisfeature in the reconstruction.

To remove the tilt, in one embodiment a line profile of thereconstructed surface along an x-axis and y-axis were compared withcorresponding profiles of the measured surface. The slopes of thereconstructed surface relative to the measured surface were estimated.Also the difference of the height of the same dominant feature (e.g.,crest) that was used for alignment of the center was determined. A planedefined by those slopes and height differences was subtracted from thereconstructed surface. In another embodiment, it has been determinedthat the tilt in the Fourier transform algorithm may come from a DCcomponent of the Fourier transform of the x and y gradients that get setto zero in the reconstruction process. Consequently, the net gradient ofthe entire wavefront is lost. Adding in a mean gradient field “untips”the reconstructed surface. As may be appreciated, such methods may beincorporated into modules of the present invention to remove the tiltfrom the reconstructions.

A comparison of reconstructed surfaces and a directly measured surfacefor a decentered +2 lens is illustrated in FIG. 6. As illustrated inFIG. 6, all of the reconstruction methods matched the surface well. TheRMS error for the reconstructions are as follows: Fourier 0.2113e−3Direct Integration 0.2912e−3 Zernike (6th order) 0.2264e−3

FIG. 7 shows a cross section of the Presbyopia Shape I reconstruction.As can be seen, the Zernike 6th order reconstruction excessively widensthe bump feature. The other reconstructions provide a better match tothe surface. The RMS error for the four reconstruction methods are asfollows: Fourier 0.1921e−3 Direct Integration 0.1849e−3 Zernike (6thorder) 0.3566e−3 Zernike (10th order) 0.3046e−3

FIG. 8 shows a cross section of Presbyopia Shape II reconstruction. Thedata is qualitatively similar to that of FIG. 7. The RMS error for thefour reconstruction methods are as follows: Fourier 0.1079e−3 DirectIntegration 0.1428e−3 Zernike (6th order) 0.1836e−3 Zernike (10th order)0.1413e−3

From the above results, it appears that the 6th order Zernikereconstructions is sufficient for smooth surfaces with features that arelarger than approximately 1-2 millimeters. For sharper features,however, such as the bump in the presbyopia shapes, the 6th orderZernike reconstruction gives a poorer match with the actual surface whencompared to the other reconstruction methods.

Sharper features or locally rapid changes in the curvature of thecorneal surface may exist in some pathological eyes and surgicallytreated eyes. Additionally, treatments with small and sharp features maybe applied to presbyopic and some highly aberrated eyes.

Applicants believe that part of the reason the Fourier transformationprovides better results is that, unlike the Zernike reconstructionalgorithms (which are defined over a circle and approximates the pupilto be a circle), the Fourier transformation algorithm (as well as thedirect integration algorithms) makes full use of the available data andallows for computations based on the actual shape of the pupil (which istypically a slight ellipse). The bandwidth of the discrete Fourieranalysis is half of the sampling frequency of the wavefront measuringinstrument. Therefore, the Fourier method may use all gradient fielddata points. Moreover, since Fourier transform algorithms inherentlyhave a frequency cutoff, the Fourier algorithms filter out (i.e., set tozero) all frequencies higher than those that can be represented by thedata sample spacing and so as to prevent artifacts from being introducedinto the reconstruction such as aliasing. Finally, because manywavefront measurement systems sample the wavefront surface on a squaregrid and the Fourier method is performed on a square grid, the Fouriermethod is well suited for the input data from the wavefront instrument.

In contrast, the Zernike methods use radial and angular terms (e.g.,polar), thus the Zernike methods weigh the central points and theperipheral points unequally. When higher order polynomials are used toreproduce small details in the wavefront, the oscillations in amplitudeas a function of radius are not uniform. In addition, for any givenpolynomial, the meridional term for meridional index value other thanzero is a sinusoidal function. The peaks and valleys introduced by thisZernike term are greater the farther one moves away from the center ofthe pupil. Moreover, it also introduces non-uniform spatial frequencysampling of the wavefront. Thus, the same polynomial term mayaccommodate much smaller variations in the wavefront at the center ofthe pupil than it can at the periphery. In order to get a good sample ofthe local variations at the pupil edge, a greater number of Zerniketerms must be used. Unfortunately, the greater number of Zernike termsmay cause over-sampling at the pupil center and introduction ofartifacts, such as aliasing. Because Fourier methods provide uniformspatial sampling, the introduction of such artifacts may be avoided.

Additional test results on clinical data are illustrated in FIGS. 9 to11. A Fourier method of reconstructing the wavefront was compared with6th order Zernike methods and a direct integration method to reconstructthe wavefront from the clinical data. The reconstructed wavefronts werethen differentiated to calculate the gradient field. The root meansquare (RMS) difference between the calculated and the measured gradientfield was used as a measure of the quality of reconstruction.

The test methods of the reconstruction were as follow: A wavefrontcorresponding to an eye with a large amount of aberration wasreconstructed using the three algorithms (e.g., Zernike, Fourier, anddirect integration). The pupil size used in the calculations was a 3 mmradius. The gradient field of the reconstructed wavefronts were comparedagainst the measured gradient field. The x and y components of thegradient at each sampling point were squared and summed together. Thesquare root of the summation provides information about the curvature ofthe surface. Such a number is equivalent to the average magnitude of thegradient multiplied by the total number of sampling points. For example,a quantity of 0 corresponds to a flat or planar wavefront. The ratio ofthe RMS deviation of the gradient field with the quantity gives ameasure of the quality of reconstruction. For example, the smaller theratio, the closer the reconstructed wavefront is to the directlymeasured wavefront. The ratio of the RMS deviations (described supra)with the quantity of the different reconstructions are as follows:Zernike (6th order) 1.09 Zernike (10th order) 0.82 Direct integration0.74 Fourier 0.67

FIG. 9 illustrates a vector plot of the difference between thecalculated and measured gradient field. The Zernike plot (noted by “Zfield”) is for a reconstruction using terms up to the 10th order. FIG.11 illustrates that the Zernike reconstruction algorithm using terms upto 6th order is unable to correctly reproduce small and sharp featureson the wavefront. As shown in FIG. 10, Zernike algorithm up to the 10thorder term is better able to reproduce the small and sharp features. Asseen by the results, the RMS deviation with the measured gradient isminimum for the Fourier method.

The mathematical development for the determination of an optical surfacemodel using an iterative Fourier transform algorithm according to oneembodiment of the present invention will now be described.

In wavefront technology, an optical path difference (OPD) of an opticalsystem such as a human eye can be measured. There are differenttechniques in wavefront sensing, and Hartmann-Shack wavefront sensinghas become a popular technique for the measurement of ocularaberrations. A Hartmann-Shack device usually divides an aperture such asa pupil into a set of sub-apertures; each corresponds to one areaprojected from the lenslet array. Because a Hartmann-Shack devicemeasures local slopes (or gradients) of each sub-aperture, it may bedesirable to use the local slope data for wavefront reconstruction.

Assuming that W(x,y) is the wavefront to be reconstructed, the localslope of W(x,y) in x-axis will be$\frac{\partial{W\left( {x,y} \right)}}{\partial x}$and in y-axis will be$\frac{\partial{W\left( {x,y} \right)}}{\partial y}.$Assuming further that c(u,v) is the Fourier transform of W(x,y), thenW(x,y) will be the inverse Fourier transform of c(u,v).Therefore, we haveW(x, y)=∫∫c(u, v)exp[i2π(ux+vy)]dudv,  (19)where c(u, v) is the expansion coefficient. Taking a partial derivativeof x and y, respectively in Equation (19), we have $\begin{matrix}\left\{ \begin{matrix}{\frac{\partial{W\left( {x,y} \right)}}{\partial x} = {i\quad 2\quad\pi{\int{\int{{{uc}\left( {u,v} \right)}{\exp\left\lbrack {{\mathbb{i}}\quad 2\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}{\mathbb{d}u}{\mathbb{d}v}}}}}} \\{\frac{\partial{W\left( {x,y} \right)}}{\partial y} = {i\quad 2\quad\pi{\int{\int{{{vc}\left( {u,v} \right)}{\exp\left\lbrack {{\mathbb{i}}\quad 2\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}{\mathbb{d}u}{\mathbb{d}v}}}}}}\end{matrix} \right. & (20)\end{matrix}$

Denoting c_(u) to be the Fourier transform of the x-derivative of W(x,y)and c_(v) to be the Fourier transform of the y-derivative of W(x,y).From the definition of Fourier transform, we have $\begin{matrix}\left\{ \begin{matrix}{{c_{u}\left( {u,v} \right)} = {\int{\int{\frac{\partial{W\left( {x,y} \right)}}{\partial x}{\exp\left\lbrack {{- {\mathbb{i}2}}\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}{\mathbb{d}x}{\mathbb{d}y}}}}} \\{{c_{v}\left( {u,v} \right)} = {\int{\int{\frac{\partial{W\left( {x,y} \right)}}{\partial y}{\exp\left\lbrack {{- {\mathbb{i}2}}\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}{\mathbb{d}x}{\mathbb{d}y}}}}}\end{matrix} \right. & (21)\end{matrix}$

Equation (21) can also be written in the inverse Fourier transformformat as $\begin{matrix}\left\{ \begin{matrix}{\frac{\partial{W\left( {x,y} \right)}}{\partial x} = {\int{\int{{c_{u}\left( {u,v} \right)}{\exp\left\lbrack {{\mathbb{i}2}\quad{\pi\left( {{ux} + {vy}} \right)}} \right)}{\mathbb{d}u}{\mathbb{d}v}}}}} \\{\frac{\partial{W\left( {x,y} \right)}}{\partial y} = {\int{\int{{c_{v}\left( {u,v} \right)}{\exp\left\lbrack {{\mathbb{i}2}\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}{\mathbb{d}u}{\mathbb{d}v}}}}}\end{matrix} \right. & (22)\end{matrix}$

Comparing Equations (20) and (22), we obtainc _(u)(u, v)=i2πuc(u, v)  (23)c _(v)(u, v)=i2πvc(u, v)  (24)

If we multiple u in both sides of Equation (23) and v in both sides ofEquation (24) and sum them together, we getuc _(u)(u, v)+vc _(v)(u, v)=i2π(u ² +v ²)c(u, v).  (25)

From Equation (25), we obtain the Fourier expansion coefficients as$\begin{matrix}{{c\left( {u,v} \right)} = {{- i}\frac{{{uc}_{u}\left( {u,v} \right)} + {{vc}_{v}\left( {u,v} \right)}}{2{\pi\left( {u^{2} + v^{2}} \right)}}}} & (26)\end{matrix}$

Therefore, the Fourier transform of the wavefront can be obtained as$\begin{matrix}{{c\left( {u,v} \right)} = {- {\frac{i}{2{\pi\left( {u^{2} + v^{2}} \right)}}\left\lbrack {{u{\int{\int{\frac{\partial{W\left( {x,y} \right)}}{\partial x}{\exp\left\lbrack {{- {i2}}\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}}}}} + {v{\int{\int{\frac{\partial{W\left( {x,y} \right)}}{\partial y}{\exp\left\lbrack {{- {i2}}\quad{\pi\left( {{ux} + {vy}} \right)}} \right\rbrack}}}}}} \right\rbrack}}} & (27)\end{matrix}$

Hence, taking an inverse Fourier transform of Equation (27), we obtainedthe wavefront asW(x, y)=∫∫c(u, v)exp[i2π(ux+vy)]dudv.  (28)

Equation (28) is the final solution for wavefront reconstruction. Thatis to say, if we know the wavefront slope data, we can calculate thecoefficients of Fourier series using Equation (27). With Equation (28),the unknown wavefront can then be reconstructed. In the Hartmann-Shackapproach, a set of local wavefront slopes is measured and, therefore,this approach lends itself to the application of Equation (27).

In some cases, however, the preceding wavefront reconstruction approachmay be limited to unbounded functions. To obtain a wavefront estimatewith boundary conditions (e.g. aperture bound), applicants havediscovered that an iterative reconstruction approach is useful. First,the above approach can be followed to provide an initial solution, whichgives function values to a square grid larger than the functionboundary. This is akin to setting the data points to a small non-zerovalue as further discussed below. The local slopes of the estimatedwavefront of the entire square grid can then be calculated. In the nextstep, all known local slope data, i.e., the measured gradients from aHartmann-Shack device, can overwrite the calculated slopes. Based on theupdated slopes, the above approach can be applied again and a newestimate of the wavefront can be obtained. This procedure is repeateduntil either a pre-defined number of iterations is reached or apredefined criterion is satisfied.

Three major algorithms have been used in implementing Fourierreconstruction in WaveScan® software. These algorithms are the basis forimplementing the entire iterative Fourier reconstruction. The firstalgorithm is the iterative Fourier reconstruction itself. The secondalgorithm is for the calculation of refraction to display in a WaveScan®device. And the third algorithm is for reporting the root-mean-square(RMS) errors.

A. Wavefront Surface Reconstruction

An exemplary iterative approach is illustrated in FIG. 12. The approachbegins with inputting optical data from the optical tissue system of aneye. Often, the optical data will be wavefront data generated by awavefront measurement device, and will be input as a measured gradientfield 200, where the measured gradient field corresponds to a set oflocal gradients within an aperture. The iterative Fourier transform willthen be applied to the optical data to determine the optical surfacemodel. This approach establishes a first combined gradient field 210,which includes the measured gradient field 200 disposed interior to afirst exterior gradient field. The first exterior gradient field cancorrespond to a plane wave, or an unbounded function, that has aconstant value W(x,y) across the plane and can be used in conjunctionwith any aperture.

In some cases, the measured gradient field 200 may contain missing,erroneous, or otherwise insufficient data. In these cases, it ispossible to disregard such data points, and only use those values of themeasured gradient field 200 that are believed to be good whenestablishing the combined gradient field 210. The points in the measuredgradient field 200 that are to be disregarded are assigned valuescorresponding to the first exterior gradient field. By applying aFourier transform, the first combined gradient field 210 is used toderive a first reconstructed wavefront 220, which is then used toprovide a first revised gradient field 230.

A second combined gradient field 240 is established, which includes themeasured gradient field 200 disposed interior to the first revisedgradient field 230. Essentially, the second exterior gradient field isthat portion of the first revised gradient field 230 that is notreplaced with the measured gradient field 200. In a manner similar tothat described above, it is possible to use only those values of themeasured gradient field 200 that are believed to be valid whenestablishing the second combined gradient field 240. By applying aFourier transform, the second combined gradient field 240 is used toderived a second reconstructed wavefront 250. The second reconstructedwavefront 250, or at least a portion thereof, can be used to provide afinal reconstructed wavefront 290. The optical surface model can then bedetermined based on the final reconstructed wavefront 290.

Optionally, the second combined gradient field can be further iterated.For example, the second reconstructed wavefront 250 can be used toprovide an (n−1)^(th) gradient field 260. Then, an (n)^(th) combinedgradient field 270 can be established, which includes the measuredgradient field 200 disposed interior to the (n−1)^(th) revised gradientfield 260. Essentially, the (n)^(th) exterior gradient field is thatportion of the (n−1)^(th) revised gradient field 260 that is notreplaced with the measured gradient field 200. By applying a Fouriertransform, the (n)^(th) combined gradient field 270 is used to derivedan (n)^(th) reconstructed wavefront 280. The (n)^(th) reconstructedwavefront 280, or at least a portion thereof, can be used to provide afinal reconstructed wavefront 290. The optical surface model can then bedetermined based on the final reconstructed wavefront 290. In practice,each iteration can bring each successive reconstructed wavefront closerto reality, particularly for the boundary or periphery of the pupil oraperture.

Suppose the Hartmann-Shack device measures the local wavefront slopesthat are represented as dZx and dZy, where dZx stands for the wavefrontslopes in x direction and dZy stands for the wavefront slopes in ydirection. In calculating the wavefront estimates, it is helpful to usetwo temporary arrays cx and cy to store the local slopes of theestimated wavefront w. It is also helpful to implement the standardfunctions, such as FFT, iFFT, FFTShift and iFFTShift.

An exemplary algorithm is described below:

-   -   1. Set a very small, but non-zero value to data points where        there is no data representation in the measurement (from        Hartmann-Shack device) (mark=1.2735916e−99)    -   2. Iterative reconstruction starts for 10 iterations        -   a. for the original data points where gradient fields not            equal to mark, copy the gradient fields dZx and dZy to the            gradient field array cx, and cy        -   b. calculate fast Fourier transform (FFT) of cx and cy,            respectively        -   c. quadrant swapping (FFTShift) of the array obtained in            step b        -   d. calculate c(u,v) according to Equation (26)        -   e. quadrant swapping (iFFTShift) of the array obtained in            step d        -   f. inverse Fourier transform (iFFT) of the array obtained in            step e        -   g. calculate updated surface estimate w (real part of the            array from step e)        -   h. calculate updated gradients from w (derivative of w to x            and y)        -   i. when the number of iterations equals to 10, finish    -   3. Calculate average gradients using the estimates from Step 2.h    -   4. Subtract the average gradients from gradient fields obtained        from Step 2.h to take off tip/tilt component    -   5. Apply Step 2.b-g to obtain the final estimate of wavefront

B. Wavefront Refraction Calculation

When the wavefront is constructed, calculation of wavefront refractionmay be more difficult than when Zernike reconstruction is used. Thereason is that once the Zernike coefficients are obtained with Zernikereconstruction, wavefront refraction can be calculated directly with thefollowing formulae: $\begin{matrix}{{C = {- \frac{4\sqrt{6}\sqrt{\left( c_{2}^{- 2} \right)^{2} + \left( c_{2}^{2} \right)^{2}}}{R^{2}}}},} & (29) \\{{S = {{- \frac{4\sqrt{3}c_{2}^{0}}{R^{2}}} - {0.5C}}},} & (30) \\{\theta = {\frac{1}{2}{{\tan^{- 1}\left\lbrack \frac{c_{2}^{2}}{c_{2}^{- 2}} \right\rbrack}.}}} & (31)\end{matrix}$where c₂ ⁻² c₂ ⁰ and c₂ ² stand for the three Zernike coefficients inthe second order, S stands for sphere, C stands for cylinder and θ forcylinder axis. However, with Fourier reconstruction, none of the Fouriercoefficients are related to classical aberrations. Hence, a Zernikedecomposition is required to obtain the Zernike coefficients in order tocalculate the refractions using Equations (29)-(31).

Zernike decomposition tries to fit a surface with a set of Zernikepolynomial functions with a least squares sense, i.e., the root meansquare (RMS) error after the fit will be minimized. In order to achievethe least squares criterion, singular value decomposition (SVD) can beused, as it is an iterative algorithm based on the least squarescriterion.

Suppose the wavefront is expressed as a Zernike expansion as$\begin{matrix}{{{W\left( {r,\theta} \right)} = {\sum\limits_{i = 0}^{N}\quad{c_{i}{Z_{i}\left( {r,\theta} \right)}}}},} & (32)\end{matrix}$or in matrix form when digitized asW=Z·c,  (33)where W is the 2-D M×M matrix of the wavefront map, Z is the M×M×Ntensor with N layers of matrix, each represents one surface of aparticular Zernike mode with unit coefficient, and c is a column vectorcontaining the set of Zernike coefficients.

Given the known W to solve for c, it is straightforward if we obtain theso-called generalized inverse matrix of Z asc=Z ⁺ ·W.  (34)

A singular value decomposition (SVD) algorithm can calculate thegeneralized inverse of any matrix in a least squares sense. Therefore,if we haveZ=U·w·V ^(T),  (35)then the final solution of the set of coefficients will bec=V·w ⁻¹ ·U ^(T) ·W.  (36)

One consideration in SVD is the determination of the cutoff eigen value.In the above equation, w is a diagonal matrix with the elements in thediagonal being the eigen values, arranged from maximum to minimum.However, in many cases, the minimum eigen value is so close to zero thatthe inverse of that value can be too large, and thus it can amplify thenoise in the input surface matrix. In order to prevent the problem ofthe matrix inversion, it may be desirable to define a condition number,r, to be the ratio of the maximum eigen value to the cutoff eigen value.Any eigen values smaller than the cutoff eigen value will not be used inthe inversion, or simply put zero as the inverse. In one embodiment, acondition number of 100 to 1000 may be used. In another embodiment, acondition number of 200 may be used.

Once the Zernike decomposition is implemented, calculation of sphere,cylinder as well as cylinder axis can be obtained using Equations(29)-(31). However, the refraction usually is given at a vertexdistance, which is different from the measurement plane. Assuming dstands for the vertex distance, it is possible to use the followingformula to calculate the new refraction (the cylinder axis will notchange): $\begin{matrix}{C^{\prime} = {\frac{\begin{matrix}{S^{\prime} = \frac{S}{1 + {dS}}} \\{S + C}\end{matrix}}{1 + {d\left( {S + C} \right)}} - {S^{\prime}.}}} & (19)\end{matrix}$

The algorithm can be described as the following:

-   -   1. Add pre-compensation of sphere and cylinder to the wavefront        estimated by iterative Fourier reconstruction algorithm    -   2. Decomposition of surface from Step 1 to obtain the first five        Zernike coefficients    -   3. Apply Equations (29)-(31) to calculate the refractions    -   4. Readjust the refraction to a vertex distance using Equation        (37)    -   5. Display the refraction according to cylinder notation

C. Wavefront Root-Mean-Square (RMS) Calculation

Finally, the wavefront root-mean-square (RMS) error can be calculated.Again, with the use of Zernike reconstruction, calculation of RMS erroris straightforward. However, with iterative Fourier reconstruction, itmay be more difficult, as discussed earlier. In this case, the Zernikedecomposition may be required to calculate the wavefront refraction andthus is available for use in calculating the RMS error.

For RMS errors, three different categories can be used: low order RMS,high order RMS as well as total RMS. For low order RMS, it is possibleto use the following formula:lo.r.m.s.={square root}{square root over (c ₃ ² +c ₄ ² +c ₅ ²)}  (38)where c₃, c₄ and c₅ are the Zernike coefficients of astigmatism,defocus, and astigmatism, respectively. For the high order RMS, it ispossible to use the entire wavefront with the formula $\begin{matrix}{{{ho}.r.m.s.} = \sqrt{\frac{\sum\limits_{n}\left( {v_{i} - \overset{\_}{v}} \right)^{2}}{n}\quad}} & (39)\end{matrix}$where v_(i) stands for the wavefront surface value at the ith location,and {overscore (v)} stands for the average wavefront surface valuewithin the pupil and n stands for the total number of locations withinthe pupil. To calculate the total RMS, the following formula may beused.r.m.s.={square root}{square root over (lo.r.m.s. ² +ho.r.m.s. ²)}  (40)

The algorithm is

-   -   1. For low order RMS, use Equation (38)    -   2. For high order RMS, use Equation (39)    -   3. For total RMS, use Equation (40)        Convergence

Convergence can be used to evaluate the number of iterations needed inan iterative Fourier transform algorithm. As noted earlier, an iterativeFourier reconstruction algorithm works for unbounded functions. However,in the embodiment described above, Equations (27) and (28) may notprovide an appropriate solution because a pupil function was used as aboundary. Yet with an iterative algorithm according to the presentinvention, it is possible to obtain a fair solution of the boundedfunction. Table 1 shows the root mean square (RMS) values afterreconstruction for some Zernike modes, each having one micron RMS input.TABLE 1 RMS value obtained from reconstructed wavefront. Real is for awavefront with combined Zernike modes with total of 1 micron error.#iteration Z3 Z4 Z5 Z6 Z7 Z10 Z12 Z17 Z24 Real 1 0.211 0.986 0.284 0.2471.772 0.236 0.969 1.995 0.938 0.828 2 0.490 0.986 0.595 0.538 1.3530.518 0.969 1.522 0.938 0.891 5 0.876 0.986 0.911 0.877 1.030 0.8610.969 1.069 0.938 0.966 10 0.967 0.986 0.956 0.943 0.987 0.935 0.9690.982 0.938 0.979 20 0.981 0.986 0.962 0.955 0.982 0.951 0.969 0.9680.938 0.981 50 0.987 0.986 0.966 0.963 0.980 0.960 0.969 0.963 0.9380.981

As an example, FIG. 13 shows the surface plots of wavefrontreconstruction of an astigmatism term (Z3) with the iterative Fouriertechnique with one, two, five, and ten iterations, respectively. For amore realistic case, FIG. 14 shows surface plots of wavefrontreconstruction of a real eye with the iterative Fourier technique withone, two, five, and ten iterations, respectively, demonstrating that itconverges faster than single asymmetric Zernike terms. Quite clearly 10iterations appear to achieve greater than 90% recovery of the input RMSerrors with Zernike input, however, 5 iterations may be sufficientunless pure cylinder is present in an eye.

Extrapolation

Iterative Fourier transform methods and systems can account for missing,erroneous, or otherwise insufficient data points. For example, in somecases, the measured gradient field 200 may contain deficient data. Isthese cases, it is possible to disregard such data points whenestablishing the combined gradient field 210, and only use those valuesof the measured gradient field 200 that are believed to be good.

A research software program called WaveTool was developed for use in thestudy. The software was written in C++ with implementation of theiterative Fourier reconstruction carefully tested and results comparedto those obtained with Matlab code. During testing, the top row, thebottom row, or both the top and bottom rows were assumed to be missingdata so that Fourier reconstruction had to estimate the gradient fieldsduring the course of wavefront reconstruction. In another case, one ofthe middle patterns was assumed missing, simulating data missing due tocorneal reflection. Reconstructed wavefronts with and withoutpre-compensation are plotted to show the change. At the same time, rootmean square (RMS) errors as well as refractions are compared. Eachwavefront was reconstructed with 10 iterations.

Only one eye was used in the computation. The original H-S patternconsists of a 15×15 array of gradient fields with a maximum of a 6 mmpupil computable. When data are missing, extrapolation is useful tocompute the wavefront for a 6 mm pupil when there are missing data.Table 2 shows the change in refraction, total RMS error as well assurface RMS error (as compared to the one with no missing data) for afew missing-data cases.

The measured gradient field can have missing edges in the verticaldirection, because CCD cameras typically are rectangular in shape.Often, all data in the horizontal direction is captured, but there maybe missing data in the vertical direction. In such cases, the measuredgradient field may have missing top rows, bottom rows, or both. TABLE 2Comparison of refraction and RMS for reconstructed wavefront withmissing data. Case Rx Total RMS RMS Error No missing data−2.33DS/−1.02DC × 170°@12.5 3.77 μm — Missing top row −2.33DS/−1.03DC ×170°@12.5 3.78 μm 0.0271 μm Missing bottom row −2.37DS/−0.97DC ×169°@12.5 3.75 μm 0.0797 μm Missing top and bottom −2.37DS/−0.99DC ×170°@12.5 3.76 μm 0.0874 μm Missing one point −2.33DS/−1.02DC ×170°@12.5 3.77 μm 0.0027 μm Missing four points −2.32DS/−1.03DC ×170°@12.5 3.76 μm 0.0074 μm

FIG. 15 shows the reconstructed wavefronts with and withoutpre-compensation for different cases of missing data. The top row showswavefront with pre-compensation and the bottom row shows wavefrontwithout pre-compensation. The following cases are illustrated: (a) Nomissing data; (b) missing top row; (c) missing bottom row; (d) missingboth top and bottom rows; (e) missing a middle point; (f) missing fourpoints. The results appear to support that missing a small amount ofdata is of no real concern and that the algorithm is able to reconstructa reasonably accurate wavefront.

With 10 iterations, the iterative Fourier reconstruction can providegreater than 90% accuracy compared to input data. This approach also canbenefit in the event of missing measurement data either inside the pupildue to corneal reflection or outside of the CCD detector.

A variety of modifications are possible within the scope of the presentinvention. For example, the Fourier-based methods of the presentinvention may be used in the aforementioned ablation monitoring systemfeedback system for real-time intrasurgical measurement of a patient'seye during and/or between each laser pulse. The Fourier-based methodsare particularly well suited for such use due to their measurementaccuracy and high speed. A variety of parameters, variables, factors,and the like can be incorporated into the exemplary method steps orsystem modules. While the specific embodiments have been described insome detail, by way of example and for clarity of understanding, avariety of adaptations, changes, and modifications will be obvious tothose of skill in the art. Although the invention has been describedwith specific reference to a wavefront system using lenslets, othersuitable wavefront systems that measure angles of light passing throughthe eye may be employed. For example, systems using the principles ofray tracing aberrometry, tscherning aberrometry, and dynamic skiascopymay be used with the current invention. The above systems are availablefrom TRACEY Technologies of Bellaire, Tex., Wavelight of Erlangen,Germany, and Nidek, Inc. of Fremont, Calif., respectively. The inventionmay also be practiced with a spatially resolved refractometer asdescribed in U.S. Pat. Nos. 6,099,125; 6,000,800; and 5,258,791, thefull disclosures of which are incorporated herein by reference.Treatments that may benefit from the invention include intraocularlenses, contact lenses, spectacles and other surgical methods inaddition to refractive laser corneal surgery. Therefore, the scope ofthe present invention is limited solely by the appended claims.

1. A method of determining an optical surface model for an opticaltissue system of an eye, the method comprising: a) inputting opticaldata from the optical tissue system of the eye, the optical datacomprising set of local gradients; b) establishing a first combinedgradient field based on the set of local gradients; c) deriving a firstreconstructed wavefront by applying a Fourier transform to the firstcombined gradient field; d) providing a first revised gradient fieldbased on the first reconstructed wavefront; e) establishing a secondcombined gradient field based on the first revised gradient field; f)deriving a second reconstructed wavefront by applying the Fouriertransform to the second combined gradient field; and g) determining theoptical surface model based on the second reconstructed wavefront. 2.The method of claim 1, wherein the optical data comprises wavefrontdata.
 3. The method of claim 2, wherein the wavefront data comprises aset of local gradients within an aperture.
 4. The method of claim 3,wherein the aperture comprises a pupil of the eye.
 5. The method ofclaim 1, wherein: the first combined gradient field comprises a firstexterior gradient field and a measured gradient field, such that themeasured gradient field is disposed interior to the first exteriorgradient field, and the measured gradient field corresponds to the setof local gradients; and the second combined gradient field comprises asecond exterior gradient field and a measured gradient field, such thatthe second exterior gradient field corresponds to the first revisedgradient field, and the measured gradient field is disposed interior tothe second exterior gradient field.
 6. The method of claim 5, furthercomprising selecting at least a portion of the second reconstructedwavefront to provide a final reconstructed wavefront, and determiningthe optical surface model based on the final reconstructed wavefront. 7.The method of claim 1, further comprising iteratively performing steps(d) through (f) to derive an n^(th) reconstructed wavefront by applyingthe Fourier transform to an n^(th) combined gradient field; selecting atleast a portion of the n^(th) reconstructed wavefront to provide thefinal reconstructed wavefront; and determining the optical surface modelbased on the final reconstructed wavefront.
 8. A method of mapping awavefront error of an eye, the method comprising: determining an opticalsurface model according to the method of claim 1; and mapping thewavefront error of the eye based on the optical surface model.
 9. Amethod of computing a correction ablation pattern for an optical tissuesystem of an eye, the method comprising: determining an optical surfacemodel according to the method of claim 1; and computing the correctionablation pattern based on the optical surface model.
 10. The method ofclaim 9, wherein: the first combined gradient field comprises a firstexterior gradient field and a measured gradient field, such that themeasured gradient field is disposed interior to the first exteriorgradient field, and the measured gradient field corresponds to the setof local gradients; and the second combined gradient field comprises asecond exterior gradient field and a measured gradient field, such thatthe second exterior gradient field corresponds to the first revisedgradient field, and the measured gradient field is disposed interior tothe second exterior gradient field.
 11. The method of claim 10, furthercomprising selecting at least a portion of the second reconstructedwavefront to provide a final reconstructed wavefront, and determiningthe optical surface model based on the final reconstructed wavefront.12. The method of claim 11, further comprising iteratively performingsteps (d) through (f) to derive an nth reconstructed wavefront byapplying the Fourier transform to an nth combined gradient field;selecting at least a portion of the nth reconstructed wavefront toprovide the final reconstructed wavefront; and determining the opticalsurface model based on the final reconstructed wavefront.
 13. A methodof modifying an optical tissue surface of an eye, the method comprising:computing a correction ablation pattern according to the method of claim9; and modifying the optical tissue surface according to the correctionablation pattern by applying a laser ablation to the eye.
 14. A systemfor determining an optical surface model for an optical tissue system ofan eye, the system comprising: a) a light source for transmitting animage through the optical tissue system; b) a sensor oriented fordetermining a set of local gradients for the optical tissue system bydetecting the transmitted image; c) a processor coupled to the sensor;and d) a memory coupled to the processor, the memory configured to storea plurality of code modules for execution by the processor, theplurality of code modules comprising: i) a module for establishing afirst combined gradient field based on the set of local gradients; ii) amodule for deriving a first reconstructed wavefront by applying aFourier transform to the first combined gradient field; iii) a modulefor providing a first revised gradient field based on the firstreconstructed wavefront; iv) a module for establishing a second combinedgradient field based on the first revised gradient field; v) a modulefor deriving a second reconstructed wavefront by applying the Fouriertransform to the second combined gradient field; and vi) a module fordetermining the optical surface model for the optical tissue system ofthe eye based on the second reconstructed wavefront.
 15. A system formapping a wavefront error of an optical tissue system of an eye, thesystem comprising the system of claim 14, wherein the plurality of codemodules further comprises a module for mapping a wavefront error of theeye based on the optical surface model.
 16. A system for computing acorrection ablation pattern for an optical tissue system of an eye, thesystem comprising the system of claim 14, wherein the plurality of codemodules further comprises a module for computing a correction ablationpattern for the eye based on the optical surface model.
 17. A system formodifying an optical tissue surface of an eye, the system comprising thesystem of claim 16, wherein the plurality of code modules furthercomprises a laser for modifying the optical tissue surface of the eyebased on the correction ablation pattern.
 18. A system for determiningan optical surface model that corresponds to an optical tissue system ofan eye, the system comprising: a) means for transmitting an imagethrough the optical tissue system; b) means, in an optical path from theimage, for determining a set of local gradients for the optical tissuesystem based on the transmitted image; c) means, coupled to the localgradient determining means, for establishing a first combined gradientfield based on the set of local gradients; d) means, coupled to thefirst combined gradient field establishing means, for deriving a firstreconstructed wavefront based on the application of a Fourier transformto the first combined gradient field; e) means, coupled to the firstreconstructed wavefront deriving means, for providing a first revisedgradient field based on the first reconstructed wavefront; f) means,coupled to the first revised gradient field providing means, forestablishing a second combined gradient field based on the first revisedgradient field; g) means, coupled to the second combined gradient fieldestablishing means, for deriving a second reconstructed wavefront basedon the application of the Fourier transform to the second combinedgradient field; and h) means, coupled to the second reconstructedwavefront deriving means, for determining the optical surface model forthe optical tissue system of the eye based on the second reconstructedwavefront.
 19. A computer program stored on a computer-readable storagemedium, the computer program comprising: a) code for transmitting animage through an optical tissue system of an eye; b) code fordetermining a set of local gradients for the optical tissue system ofthe eye based on the transmitted image; c) code for establishing a firstcombined gradient field based on the set of local gradients; d) code forderiving a first reconstructed wavefront based on the application of aFourier transform to the first combined gradient field; e) code forproviding a first revised gradient field based on the firstreconstructed wavefront; f) code for establishing a second combinedgradient field based on the first revised gradient field; g) code forderiving a second reconstructed wavefront based on the application ofthe Fourier transform to the second combined gradient field; and h) codefor determining an optical surface model for the optical tissue systemof the eye based on the second reconstructed wavefront.
 20. The computerprogram of claim 19, further comprising: code for computing a correctionablation pattern based on the optical surface model.
 21. The computerprogram of claim 20, further comprising: code for delivering a laserenergy to the eye based on the correction ablation pattern.